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ABSTRACT 


i  This  report  examines  the  Fourier  transform  method  of  restorinq 
degraded  images  of  point  objects.  The  principal  conclusion  supported 
by  this  study  is  that  the  major  problem  in  restoring  these*  imaqes  is 
the  presence  of  poles  iri  the  restored  image  spectrum.  A  series  of 

thr®e  re^tored  imr‘9e  quality  criteria  graphed  as 
functions  of  image  degradation,  for  a  computer  modeled  degraded 
imaging  system,  are  presented  to  support  this  conclusion.  Also 
spectrum  ^  3  th°r0U° '  mathemati cal  analysis  of  the  restored  image 
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CHAPTER  I 
INTRODUCTION 


A.  Statement  of  the  Problem 

We  know  that  atmospheric  turbulence  degrades  the  images  of  objects 
that  are  recorded  by  optical  systems.  Mirages  and  star  twinkling  are 
typical  examples.  It  is  important  for  some  applications  that  we  be 
able  to  extract  as  much  information  as  possible  from  these  degraded 
images  concerning  the  actual  object  especially  in  cases  where  this 
object  is  an  unknown  quantity.  Because  the  human  visual  system  is  an 
inefficient  processor  of  this  data,  it  is  necessary  to  resort  to  otner 
processing  techniques.  Out  of  tin's  need  for  more  efficient  processing 
techniques  emerqes  the  study  of  image  restoration. 

The  goal  of  this  study  is  to  specify  the  limits  over  which  one  such 
processing  technique,  the  Fourier  transform  method,  does  restore  a 
urbulence  degraded  image  so  that  it  is  a  more  accurate  representation 
o  the  actual  object.  For  this  study,  we  will  model  a  turbulence 
degraded  imaging  system,  and  apply  the  Fourier  transform  processing 

technique  to  the  resulting  degraded  images  to  determine  the  actual 
object. 

B-  Background  Discussion  of  Imaging  in  the  Presence  of 
Atmospheric  Turbulence 


Objects  which  are  imaged  through  the  atmosphere  are  degraded  as  a 
result  of  random  temperature  fluctuations  in  that  atmosphere.  In  this 
section  a  brief  discussion  of  there  degrading  effects  is  presented. 


Imaging  through  a  turbulent  atmosphere  is  not  a  new  phenomena, 
for  man  s  own  imaging  system,  the  eye,  has  always  had  to  contend  with 
this  problem  whenever  a  temperature  gradient  was  present  along  his  line 
or  signt.  One  of  the  most  familiar  instances  today  of  the  eye  imaging 
through  a  tuioulent  atmosphere  occurs  while  driving  along  a  highway 
which  is  reradiating  solar  energy.  As  one  looks  at  a  car  or  the 
scenery  some  distance  ahead,  the  objects  seem  to  be  dancing  and  are 
blurred  because  of  the  atmospheric  effects.  In  this  case  the  phenomena 
is  commonly  referred  to  by  the  phrase  "the  heat  is  rising". 

The  atmosphere  is  a  gas  and  thus  as  it  is  warmed,  its  density,  and 
•?he+ltS  in?ex  of  re  fraction,  changes.  The  warmer  layers  rise  and  mix 
with  the  cooler  and  more  dense  layers  which  are  falling  due  to  gravi¬ 
tational  forces.  This  mixing  causes  additional  refractive  index 
changes  to  occur  These  changes  are  referred  to  as  refractive  index 

obviously  random  functions  of  time  which  are  best 
described  by  statistical  means.  (See  Tatarski  [l,2j.) 
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A  light  wave  propagating  through  a  non-turbulent  atmosphere  is 
focused  by  a  telescope  to  a  specific  point  in  the  image  plane.  How¬ 
ever,  when  the  propagation  path  contains  refractive  index  fluctuations 
one  may  consider  effects  occuring  on  two  spatial  scales  which  degrade 
the  ideal  imaging  situation.  First,  for  those  refractive  index  fluc¬ 
tuations  whii.n  vary  slowly  over  a  spatial  scale  larger  than  the  tele¬ 
scope  input  aperture,  the  image  plane  effect  is  similar  to  that  which 
occurs  when  a  prism  is  placed  in  front  of  the  telescope  input  aperture. 
That  is,  the  incoming  rays  are  bent  causing  the  focused  image  to  be 
shifted  to  a  new  inage  plane  coordinate  location.  Second,  for  those 
refractive  index  fluctuations  which  vary  raoidly  over  a  spatial  scale 
comparable  to  the  input  aperture,  the  result  observed  in  the  optical 
system  image  plane  is  that  of  image  break-up. 

Both  of  these  effects  obviously  perturb  the  incoming  wave  phase 
fronts.  By  including  a  random  time  variation  in  these  phase  front 
perturbations,  the  resulting  image  plane  pattern  is  a  "dancing",  blurred 
and  broken-up  corruption  of  the  ideal  image. 

C.  Mathematical  Concepts  and  the  Fourier  Transform  Method  of 
Image  Restoration 

T^e  familiar  one-dimensional  mathematical  concepts  of  the  convo¬ 
lution  integral  and  the  Fourier  transform  pair  are  extended  to  a  two- 
dimensional  generalization.  These  two-dimensional  generalizations 
form  the  mathematical  foundation  on  which  the  Fourier  transform  method 
of  image  restoration  is  based. 

First,  consider  the  extension  of  one-dimensional  linear  system 
theory,  conmonly  associated  with  electrical  engineering,  to  the  two- 
dimensional  linear  imaging  system.  It  has  been  shown  (s?e  Goodnan  [3]) 
that  the  image  plane  intensity  function  and  the  object  plane  intensity 
function  for  an  incoherent  object  are  related  by  a  convolution  integral 
of  the  form: 


(i)  ri (x2 *y2) =  ij  h(vxi’y2-yi)  Mxi>yj)  dxi  d^i  » 


where  Ii(x2,y2)  and  I p ( x  1  »y  1 )  are  the  image  and  object  plane  intensity 
functions,  respectively.  This  integral  should  be  recognized  as  the 
extension  to  two  dimensions  of  the  one-dimensional  integral  familiar  to 
all  electrical  engineers,  which  relates  the  output  ano  input  of  a  linear 
system.  As  in  linear  network  theory,  the  two-c'mensional  generalization 
has  associated  with  it  a  unit  impulse  response  h(x2,y2).  This  impulse 
response  is  obtained  by  placing  an  intensity  point  source  in  the  object 
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PLANE 


Fig.  1.  Imaging  coordinate  systems. 

plane  This  impulse  rerpon  e  is  referred  to  as  the  intensity  point 
spread  function  (PSF)  arid  is  defined  as:  P 


(2) 


oo 

Ii(x2,y2)  =  h(x2,y2)  =  jj  h(x2-x1,y2-y1)  6(xx ,y x )  dXj  dy, 

-  CO 


on . 


where  <5^^)  is  a  two-dimensional  Dirac  delta  functi 

The  integral  in  Eq.  (1)  is  valid  only  over  that  portion  of  the 
object  plane  for  which  h(x2,y2)  is  independent  of  the  location  of 

in1  •  Z  region  over  wnich  The  integral  is  said  to  be  spatial  1 

.n van  ant  is  known  as  an  isopl  anatic  patch  and  is  analogous  to  ?he  ti¬ 
ll. variance  condition  in  electrical  circuits. 

pair  3  two'di^ mens 1 0,191  generalization  of  the  Fourier  transfm 


(3)  F(x,y)  F(KX,K^)  , 
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where  the  symbol  <-*  indicates  that  the  two  functions  are  Fourier  trans¬ 
form  mates,  we  have: 


,  ,  n  -J(Kxx  +  Kvy) 

(4)  F(KX,K  )  =  f(x,y)  e  y  dx  dy 


and 

(5)  f(x,y)  = 


1  u  r/„  ,  j(Kxx+Kyy) 

F(Kx,Ky)  e  y  dKx  dKy 


(2*)2  JJ 


In  Eqs.  (4)  and  (5),  Kx  and  Ky  are  defined  as  spatial  frequencies  with 
the  dimensions  of  radians/mm.  Applying  this  generalized  Fourier  trans¬ 
form  to  Eq.  (1),  we  see  from  the  convolution  theorem  that: 


(6a)  I.  (Kx,Ky)  =  H(Kx,Ky)  IQ(Kx,Ky) 


The  function  H(Kx,Ky)  is  known  as  the  system  modulation  transfer 
function  (MTF)  and  is  analogou'  to  the  network  transfer  function  used 
in  linear  circuit  theory.  The  functions  I* (KX,KV)  and  I0(KX,KV)  are 
the  spatial  spectra  of  the  image  and  object  plane  intensity  functions, 
respectively. 


The  recorded  image  plane  intensity  will  be  a  degraded  form  of 
Ip(xi.yi),  the  object  intensity  function,  if  the  region  between  the 
object  plane  and  the  telescops  input  aperture  is  assumed  to  be  a  turbu¬ 
lent  medium.  However,  we  can  obtain  the  MTF  of  the  system,  including 
the  turbulence  effects,  by  placing  a  point  source  in  the  same  i sop! ana- 
tic  patch  occupied  by  Io(xi,yj).  The  input  intensity  function  Io(xi,yi) 
can  then  be  determined  by  dividing  the  degraded  image  plane  intensity 
function  spectrum  by  the  system  MTF  (the  quotient 


(6b) 


I0(Kx,Ky) 


Ij  (Kx,Ky) 
H(Kx,Ky) 


is  defined  as  the  restored  image  spectrum  (RIS) )  and  taking  the  inverse 
Fourier  Transform  of  this  RIS.  The  result  of  these  operations  is  the 
restored  image,  which  is: 


<7>  ■Rtxi.yi)  *  V*i.yi)  -  nVlLIci(K*,Ky) e 

1  ,(■  (Kx,Ky)  j(Kxx+  Ky y) 


j(Kxx  +  Kyy; 


dKx  dKy 


(2")2  n  H(Kx,Ky) 

where  ^(xi.y^)  is  defined  as  the  restored  image. 


dK 


x  dKy 


wh^hFwrca/?-d°!!ly  V^in9  turbulent  medium,  an  isoplanatic  region, 
which  was  defined  earlier  as  the  region  over  which  Eq.  (1)  is  spatially 

lTvarvi™  JraHc°-t^finedflby+hharp  boundaries  but  rather  by  mom  smooth¬ 
ly  varying  transitions.  As  the  separation  between  the  object  and  the 

point  source  increases,  the  phase  front  perturbations  across  the  tele- 

srope  input  aperture  become  different  since  the  beams  have  propagated 

along  path-  haying  slightly  different  small  and  large  scale  refracti ve 

index  fluctuations.  Thus ,  one  would  expect  that  as  the  reference  point 

source  is  moved  through  such  transition  regions,  the  quality  of  the 

WOil  d  f  grease  t0  the  point  where  it  cannot  be  consider¬ 
ed  as  a  valid  restored  image. 

D .  Summary  of  ^nage  Restoration  Methods 

In  the  past,  there  have  been  various  methods  suggested  for  the 
restoration  of  turbulence  degraded  images.  Among  these  schemes  have 
been  wavefront  reconstruction  (Gaskil]  [4],  Goodman  [5]);  maximum 

°ratl0n-(F!i1^den  E6]);  sPa«al  filtering  (Mueller  [7] , 
Reynolds  [7]);  constrained  deconvolution  (MacAdams  [8]);  and  the 

Fourier  transform  method  (Harris  [9],  McGlamery  flOj)  1  These  and  other 

Effects  Stherbthan^t  image  de9radation  ^suiting  from 

effects  other  than  atmospheric  turbulence,  such  as  image  motion  or 

optical  and  electrical  degradation  (Huang  [11],  Andrews  [12], 


ac  ifTanni?IeS+nt^tUdy  1S  c°,?cerned  W1’ th  the  Fourier  transform  method 
as  it  applies  to  the  restoration  of  turbulence-degraded  imaqes  It 

sb°uld  1 \ POin+,ted  °Ut’  however’  that  thl‘s  rethod  also  applicable  to 
some  of  the  other  sources  of  image  degradation  such  as  image  motion. 

crhom!hh  pi?vi0us  w°r*5  ™volvin9  the  Fourier  transform  restoration 

tatlV6,in  natUre-  Tt  has  demonstrated  this  method 
using  artificial ly  generated  or  computer  modelled  image  degrading 

fects.  In  this  study  we  emphasize  the  quantitative  aspects  of  the 
restored  image  quality  obtained  using  the  Fourier  transform  method. 

InLi^cr  qu®1Ity-15  ,pre^ented  as  a  function  of  deviations  from  the 
ideal  PSF  which  yielas  the  perfectly  restored  image. 

The  question  we  attenpt  to  answer  in  this  study  is,  how  larqe 

t^nmHf.ro^rH1-005  before  the  Fourier  transform  method  fails 

to  produce  good  image  restoration  results?  The  results  obtained  from 

answering  this  question  may  serve  as  the  input  data  needed  to  answer 
a  *"*  What  18  the  is^a"atl’c  Of 

F .  Organization  of  the  Study 

.  .  Sn  Chapter  I,  we  have  discussed  the  image  degrading  effects  of  a 
turbu.ent  atmosphere.  We  have  described  a  typical  imaging  system  and 
presented  the  necessary  mathematical  background  required  to  apply  the 
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Fourier  transform  method  of  image  restoration  to  a  turbulence  deqraded 
image.  3 


In  Chapter  II,  we  generate  a  mathematical  model  that  describes 
the  turbulence  degraded  wavefronts  across  the  input  aperture  of  the 
imaging  system  as  a  polynomial  series.  Using  this  model  we  then 
derive  an  equation  for  the  restored  image  spectrum  (RIS).  The  actual 
restored  image  obtained  by  using  the  Fourier  transform  method  can  then 
be  determined  by  calculating  the  inverse  Fourier  transform  of  this  RIS. 
The  discrete  form  of  this  RIS  is  also  presented  in  this  chapter.  It 
is  thTs  discrete  equation  that  is  programmed  on  a  digital  computer 
which  calculates  the  inverse  Fourier  transform  (yieldinq  the  restored 
image)  for  various  model  parameters. 


In  Chapter  III  we  first  define  what  is  meant  by  the  phrase  "an 
ideally  restored  image".  Based  on  this  definition  we  then  note  that 
the  presence  of  poles  in  the  RIS  may  cause  the  inverse  Fourier  trans¬ 
form  of  the  RIS  to  be  a  very  poorly  restored  image.  As  a  result,  we 
then  conduct  an  excensive  analysis  of  the  RIS  and  its  poles,  the  out¬ 
come  being  a  series  of  equations  relating  the  number  of  poles  and  their 
positions  to  the  degraded  wavefront  model  parameters.  These  equations 
serve  as  input  data  for  a  quantitative  analysis  of  restored  imaqe 
quality  versus  RIS  pole  position  and  number.  Finally,  in  Chapter  III 
we  define  the  parameters  that  we  will  use  to  describe  the  quality  of 
a  restored  image. 


In  Chapter  IV,  quantitative  results  of  restored  image  quality 
are  presented  as  a  function  of  the  degraded  wavefront  model  parameters. 

ree  case:;  are  thoroughly  examined  that  relate  these  model  parameters 
to  specifically  known  RIS  pole  placement  and  number;  the  result  is 
the  definition  of  regions  of  successful  image  restoration  separated  by 
regions  characterized  by  gross  departures  from  the  restored  imaqe 
quality  criterion.  A  series  of  graphs  is  included  for  each  case 
showing  the  regions  of  successful  and  unsuccessful  image  restoration. 


.  chaPter  V  is  devoted  to  a  summary  of  the  entire  study.  Also,  a 
series  of  conclusions  are  drawn  relating  the  results  of  this  study  to 
the  task  of  restoring  turbulence  degraded  images 


F.  Summary 


The  goal  of  this  study,  as  stated  previously,  is  to  determine  the 
limits  over  which  the  Fourier  transform  method  of  image  restoration 
does  indeed  produce  a  valid  representation  of  the  actual  input  object 
A  secondary  goal  is  to  specify  these  limits  in  such  a  manner  that  they 

aLln?Ut  ?ata  f0r  a  future  stud*  bating  the  restored  image 
quality  to  the  isoplanatic  patch  size.  y 
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CHAPTER  II 

CALCULATION  OF  THE  RESTORED  IMAGE  SPECTRUM 


A.  Introduction 

In  Chapter  I  we  briefly  described  the  degrading  effects  of  the 
large  and  small  scale  index  fluctuations  on  an  imaging  system.  We  also 
indicated  that  the  Fourier  transform  method  will  be  the  approach  used 
in  this  study  to  restore  such  turbulence  degraded  images. 

In  this  chapter,  we  turn  our  attention  to  a  more  detailed  mathe¬ 
matical  model  of  the  imaging  system  from  which  ..he  restored  image 
spectrum  (RIS)  is  calculated  in  terms  of  phase  front  parameters.  The 
incoming  turbulence  degraded  optical  phase  front  will  be  expanded  in 
terms  of  a  general  second  degree  quadratic  equation.  Using  this  ex¬ 
pansion  we  then  calculate  the  restored  image  spectrum  for  both  con¬ 
tinuous  and  discrete  values  of  the  spatial  frequencies  and  Ky. 

These  results  will  be  used  in  Chapter  III  where  we  thoroughly  examine 
the  RIS. 


B.  The  Imaging  System  Model 

This  imaging  system  model  discussion  is  divided  into  two  sections. 
First,  the  physical  layout  and  coordinate  system  used  are  described. 
Second,  the  expansion  for  the  turbulence  degraded  input  electric  fields 
(E  fields)  is  derived. 


We  choose  the  coordinate  system  in  Fig.  2  for  the  imaging  system. 

The  (xo.yo)  plane  contains  the  telescope  input  aperture  which  is 
assumed  to  be  square,  with  width  W  =  2x  =  2y  .  This  assumption  sim¬ 
plifies  the  mathematics  but  does  not  detract  from  the  fundamental 
results . 


The  (x^y^  object  plane  is  located  a  distance  do  from  the  input 
aperture  and  contains  the  reference  point  source  and  the  object  that 
we  are  imaging.  The  medium  separating  these  two  planes  is  assumed  to 
be  a  turbulent  atmosphere  having  randomly  varying  refractive  index 
fluctuations.  The  (X2,y2)  plane  is  located  a  distance  d^  from  the 
input  aperture  and  is  the  image  plane  of  the  telescope.  The  distances 
do  and  d.j  are  conjugate  distances  since  we  are  considering  an  imaging 
situation.  The  region  between  the  (x0,y0)  and  (x2,y2)  planes  is  as¬ 
sumed  to  be  void  of  any  turbulent  effects. 

Both  the  reference  source  and  the  object  are  assumed  to  be  mono¬ 
chromatic  point  sources  which  are  described  mathematically  by  the  two- 
dimensional  Dirac  delta  function.  The  separation  between  these  two 
points  is  the  distance  ds .  The  imaging  system  is  assumed  to  be  linear. 
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Fig.  2.  Coordinate  system  for  imaging  system  model. 


thus,  analysis  of  an  extended  object  would  require  integration  of  the 
point  object  results  over  all  points  in  the  extended  object. 


The  field  distribution  over  the  input  aperture,  for  a  point  object 
located  at  a  finite  distance  d0  and  in  the  absence  of  turbulence,  is 
merely  a  spherical  wave  having  uniform  phase.  However,  in  the  presence 
of  turbulence  and  assuming  negligible  amplitude  effects,  we  use 
Frieds  [14]  assumption  that  the  perturbed  E-field  phase  fronts  may  be 
described  by  a  general  second  degree  quadratic  equation  involving  Xr 
andyQ.  That  is,  the  instantaneous  phase  $  may  be  expanded  in  the  form: 


(8a)  <j>  =  axQ2  +  bx0yQ  +  cyQ2  +  dxQ  +  eyQ  +  f 


This  formulation  can  represent  any  wavefront  that  involves  spheric?! 
and  hyperbolic  perturbations,  wavefront  tilt  and  an  average  phase  term. 

An  alternative  representation  for  the  instantaneous  phase  *  which 
clearly  defines  these  distortions  is  the  polynomial  series  used  by 
Collins  [15]. 


(8b)  4>  =  X  ap  Fn(x0,y0) 

n=l 

where: 

Fi(x0,y0)  =  i- 
w 

W2 

+  y02  -  j-)  /  (m/2)2 
-  y02)  /  (m/2)2 
^6  (xo  0  ^  =  JJ*  x<>y</  (W/2)2 


Each  polynomial  has  a  particular  interpretation.  Fj(x0,y0)  represents 
a  constant  phase  contribution;  F2(x0,y0)  and  F3(x0,y0)  represent  tilts 
in  the  xp  and  yQ  directions  respectively;  F4(x0,y0)  represents  a  spher¬ 
ical  distortion;  and  F5(xg,y0)  and  F6(x0,y0)  represent  hyperbolic 
distortions.  These  functions  are  chosen  to  be  orthonormal  over  a 


square  of  width  W.  The  coefficients  a,  bf  c,  d, 

a 


are  related  to  the  coefficients  a!,  a!,  a 


following  equations  in  which 


1’  2’  3’ 


and 

and 


f  in  Eq.  (8a) 


a'  by  the 

b 


r  _  1 3  pr  1 

Lw  '  W  7  J2  WW  * 

.  _  1  3  [T  1  .  ,  , . 

W  2  >|2  (W/2)2  u4  a 5 ' 

.  =  I _ L_  « 

"  W  (W/2)2  85 

C  =  i  1  IF _ 1 

W  2  \|2  (W/ 


(W/2)' 


W/2  W 


Cu«  +  as> 


(a4  ~  a5)  -  C^(a4  -  a5) 
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e  =  —  —  a' 
W/2  W  a3 

f  =  <ll  _  [E 
W  J2  w 


Se  Sar?nLntraof°;h?lV!.n  }"  ^  (8a)  S1!'  be  used  ““lively  throughout 
tne  reminder  of  this  study  because  of  its  computational  simplicity 

coefficients  tTfacil  IV  ar®.g1v®n  in  terms  of  ^oth  sets  of 

the  [  !  pJnlnc  i U  t?,interp^tati0n  ln  light  of  other  vrork  using 

tne  tq.  (8b)  expansion.  Thus,  using  the  expansion  in  Eq  (8a)  the  tur- 

bulence  degraded  field  distribution  over  the  aperture  plane  Resulting 
from  a  point  object  may  be  expressed  as:  resulting 


(9) 


Ea(x0,yfl)  =  e 


J  K*  +  bxQy0 


CV  + 


dx„  +  ey 
0  Jo 


f) 


s™tem'Io*?.di5tribUt1°n  15  the  genera'  fo™  assumed  for  the  Imaging 

fhaf  the  field  of  0Ptics,  it  is  a  well  known  fact  (see  Goodman  h  I) 
that  when  a  lens  system  images  an  object,  the  field  distribution  over 
the  input  aperture  and  the  field  distribution  over  the  imiqe  pi ane 

trantfom  9  ^stant  term  and  a  phase  term  from  being  exact  Fourier 
transform  mates.  This  relationship  is  expressed  as: 

00 

(10)  E.-(x„vj  .  -AA-KJcff  E  (x  v  )ej(k/di)(x°VV2>  . 

1  2’2  (2^\d.|e  Jj  ta^x0,y0)e  dx0  dyD 

*  CD 

di  strl butions °respe c ti vefy^2 ’  ^  **  aperture  and  1nH9e  plane  field 
If  we  define  the  spatial  frequencies  as: 

A  kx. 


(11a) 

and 


K  =  0 

x  ar- 


(11b)  K  £  'll 
di 


then  the  more  familiar  Fourier  transform  expression  is  obtained. 
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i  1  fr  J  KXX.«K  y. 

(12)  E^  (x?,y?)  =  A  e  ^yy  J  Ea(Kx,Ky)o  *  dKx  dKy  . 

-  u* 

By  substituting  the  above  definition  of  spatial  frequencies  into 
Eq.  (9),  the  aperture  plane  E-field  distribution  becomes: 

j  £( d-j / k ) 2  ( aKx 2+bKx Ky+cKy 2 )  +  (d-j/k )  (dKx+eKy)  +  t] 

(13)  Ea(Kx>Ky)  =  e 


In  this  study,  where  we  are  interested  in  the  instantaneous  value  of  the 
E-field  in  Eq.  (13),  we  choose  not  to  include  the  wavefront  tilt  terms 
(the  d  and  e  terms)  and  the  average  phase  term  (the  f  term)  in  the 
analysis  of  the  point  object.  The  wavefront  tilt  terms  merely  shift 
the  image  to  a  different  coordinate  location  but,  for  short  exposures, 
do  not  actually  degrade  the  image  detail  (recall  the  time  shifting 
property  in  one-dimensional  Fourier  transform  theory).  The  constant 
phase  term  is  merely  a  multiplicative  factor  that  drops  out  anyway  when 
the  intensities  are  calculated.  Thus,  the  actual  field  distribution 
assumed  for  the  imaging  system  model  is: 


(14) 


j  (( d^ /k ) 2 (aKx2  +  bKxKy  +  cKy  )J 
Ea(Kx>Ky)  =  e 


The  parameters  a,  b  and  c  are  defined  as: 

(15a)  a  =  2Trai/(W/2)2  =  2Trai/xm2 

(15b)  b  =  2TTb1/(W/2)2  =  ZTTbj/X^ 

(15c)  c  =  2-rrc  j/ (W/2)2  =  2irCl/ym2 


where  al5  bj  and  Ci  specify  the  amount  of  wavefront  distortion,  in 
fractions  of  a  wavefront,  at  the  edge  of  the  input  aperture.  Further¬ 
more,  since  the  input  aperture  is  now  considered  a  spatial  frequency 
plane,  the  aperture  dimensions  define  the  system  E-field  spatial 
frequency  band  limit  Km  to  be: 

(16)  Km  =  xm  =  Hj"  ^m  =  <£j" ^ 

Rementer  that  the  aperture  has  been  assumed  to  be  square  with  width 
^  -  ^xm  ~  ^m* 
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C •  ihe  Continuous  Restored  Image  Spectrum 

Using  the  model  discussed  in  section  B  ana  the  degraded  aperture 
plane  L-field  given  in  Eq.  (14),  the  restored  image  spectrum  (RIS)  is 
next  calculated.  ' 

The  first  step  toward  finding  the  RIS  is  the  calculation  of  the 
image  plane  intensity  function  and  its  spectrum  in  terms  of  the  E-field 
given  in  Eq.  (14).  This  image  plane  intensity  function  is 


{U)  I (x2,y2)  =  Ei(x2,y2)  E.*(x;  ,y2) 


where  the  notation  indicates  the  complex  conjugate  cf  the  function. 
(Note  that  Zq,  the  impedance  of  free  space,  has  been  dropped.  It  is 
merely  a  multiplicative  constant  that  would  be  cancelled  in  later 
steps  where  ratios  of  intensities  are  forred.)  E^x^y  )  is  given  by: 


(18)  Ei(x2 ,y2)  =  jjw(Kx,Ky)  Ea(Kx,Ky)  e(  X  2  ^  dKx  dKy 

-oo 

where  W(KX,|C)  is  the  aperture  function  which  takes  into  account  the 
finite  sizecr  input  aperture  The  function  W(Kx,Ky)  is  defined  as: 


1 


~£m  -  jjx1  Km 

‘Km1  Ky<  Kjh 


(19)  W(Kx,Ky)  =  - 

-0  otherwise 

From  Fourier  transform  theory  we  know  that  if: 
Ei (xf ,y2)  <->  Ea(Kx,Ky) 


then 


Ei*(x2,y2)  Ea*(-Kx«-Ky) 


and,  additionally,  that  multiplication  in  the  spatial  plane  is  equivalenl 
to  convolution  in  the  spatial  frequency  plane.  Using  these  facts,  the 
image  plane  intensity  is  found  to  be  the  Fourier  transform  of  the  con¬ 
volution  integral  given  in  Eq.  (20)  below. 


12 


(20)  I(x2,y2)  =  Ei(x2,y2)  Ei*(x2,y2) 


1  '  °° 

T27F /[  w(Kx-kx*Ky"ky)  w*^kx*-ky) 

00 

^a(kx"kx»ky“ky)  ^a^“kx»“ky )  ^kx  ^ky 

The  variables  kx  and  k  are  dummy  variables  that  have  the  dimensions  of 
spatial  frequency.  3 

By  letting  kx  =  -kx  and  ky  =  -ky  and  multiplying  Eq.  (20)  by  (2t)2, 
we  define  the  image  plane  function  lj(x2,y?)  as: 

(2D  !i  (x2 >y2)  =  (2tt)2I  (x2  ,y2)  <->  I-j  (KX,K  )  =  C(KX,KJ  = 


/[  W(Kx+kx,Ky+ky)  W*(kx,ky) 


Ea(Kx+kx,Ky+ky)  Ea*(kx,ky)  dkx  dky 

(kx’U  )l  the  sPatial  frequency  spectrum  of  Mx-.yJ.  The 
in  Eq.  (21)  has  the  form  of  a  correlation  function  and  hence 
is  referred  to  as  the  aperture  plane  correlation  function,  which  is 
defined  as  C(Kx,Ky). 

Substituting  Eq.  (14)  into  Eq.  (21),  we  have: 


where  Ij 
integral 


(22)  I-j  (x?  ,y?)  C(Kx,Ky)  =  Jj  W(Kx+kx,Ky+ky)  W*(kx,ky) 

"00 

ej(<f,-A)2[a(Vkx)2  +  b(Kx+kx)(Ky+ky)  +c(Ky+ky)2] 

j(d^/k)2[ak  2  +  bkxky  +  ckv2] 

e  y  y  dkv  dk 

x  y 

When  the  phase  term  involving  only  k  and  Kv  is  brought  outside  the 
integral,  we  have:  y 
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j(dj/k)2  (aKx  2+  bKxKy  +  cKy2  ) 


(23)  Ii  (x?,y2)  C(Kx,Ky)  =  e 

II  w(Kxnx,Vky;  w*(kx,ky)  ej(dl/k)? 


dkx  dky  . 


The  details  of  evaluating  the  above  integral  are  carried  out  in 
Appendix  A.  T^e  final  result  of  this  calculation  is  presented  below 


(24)  Ij(x2,y?)  C(KX,KV)  = 


(2Km  -  IM) 


(2Km  -  IM) 


Sinf  (di/k)2( 2aKx+bKyh( 2Km- 1 Kx | ) } 


/k)2(2aKx+bKy)^(2K  -|KJ) 


TSinf  ( 

Lldi 

[Sinf  ( d-j/k)2 (2cKy+bKx)55(2Kfn- 1 Ky  | ) } 
L  (di/k)2(2cK  y+bKx)*s(2Km-|Kv|) 


Note  that  C(Kx,Ky),  being  tho  transform  of  an  image  intensity  function, 
has  a  spatial  frequency  bandwidth  twice  that  of  the  image  plane  E-field 
given  in  Eq.  (16).  The  region  over  which  C(K  ,K  )  is  defined  is: 


~2Km  -  Kx 


2Km 


-2Km  <  Ky<  '2Km 

(24)j  nthe  Pr°duct  of  two  triangular  functions  (2Km  -  |KJ)  and 
).2Km  “  j^y N  each  of  which  is  modulated  by  a  function  having  a  Sin  x/x 
format  (ctmmonly .referred  to  as  a  Sine  function).  The  degrading  effects 
of  the  atmospheric  turbulence,  as  represented  by  the  parameters  a,  b, 
and  c,  are  introduced  through  the  arguments  of  the  Sine  functions.’ 

Note  that  when  a  =  b  =  c  =  0  (i  .e. ,  in  the  absence  of  turbulence),  the 
correlation  function  is: 


(25)  C(Kx,Ky)  =  (2K,,,  -  | Kx | )  (21^  -  |Ky|) 

This  function  is  the  two-dimensional  generalization  of  the  familiar 
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correlation  of  two  one-dimensional  rectangular  pulses.  For  the  two- 
dimensional  case,  Cg.  (25)  is  the  correlation  of  the  aperture  funcrion 
W(Kx,Ky)  with  itself. 

The  reference  source  is  also  a  point  source;  thus,  Eq.  (14)  also 
describes  the  resulting  aperture  phase  distribution  when  there  is  zero 
separation  (d_  =  0)  between  the  reference  point  and  object  point. 
However,  for  finite  separations  between  these  two  points,  the  atmos¬ 
pheric  effects  on  the  waves  radiating  from  the  two  points  will,  in 
general,  be  different.  These  differences  are  incorporated  in  thr 
aperture  plane  expansion  for  the  reference  E-field  by  defining  the 
reference  phase  parameters  as: 


(26a)  a0  =  a  +  Aa 
(26b)  b0  =  b  +  Ab 
(26c)  c0  =  c  +  ac 


where  Aa,  Ab,  and  ac  represent  small  perturbations  in  the  parameters 
a,  b,  and  c,  which  describe  the  aperture  plane  E-field  for  the  point 
object. 


The  reference  point  aperture  plane  field  distribution  is  obtained 
by  substituting  Eqs.  (26a),  (26b),  and  (26c)  into  Eq.  (14)  to  obtain: 


(27)  Ea(Kx,Ky)  =e 


=  e 


j(d1-/k)2(a0Kx2  +  boKxKy  +  cQKy2) 

j( d1  / k f  ((a+Aa)Kx2  +  (b+Ab)KxK  +  (c+ac)K  2) 


The  modulation  transfer  function  (MTF)  is  the  system  transfer 
function  obtained  when  the  reference  input  is  a  point  source  (see 
Chapter  I).  Actually  we  have  already  derived  this  function  in  Eq.  (24) 
where  the  object  was  also  assumed  to  be  a  point  source.  Thus  the  MTF 
is  obtained  by  substituting  Eqs.  (26a),  (26b),  and  (26c)  into  Eq.  (24). 
Remember  that  the  resulting  function  is  the  Fourie*"  transform  of  the 
system  PSF  (see  Chapter  I).  The  MTF  is: 
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(28) 


H(VKy)  =  (2Km-lKxl> 


Sinf  (di/k)y(2  |(a+Aa)Kx+(b*Ab )Ky[ '^(2Km- [ Ky  | ) ) 
.  (d-f/k )? (2  [(a+Aa)Kx+(b+Ab)KyJJs(2Km-|Kx|) 


SH_(  V k ) ■ 2  ( 2  [ ( c+ A c ) Ky+  ( b+A b ) K J  h ( 2 Km- 1  Ky  | ) } 
.  (d1-/k)2(2[(c+Ac)Kyx(b+Ab)Kx]js(2Km-|Ky!) 


We  now  have  calculated  the  spectrum  for  the  image  plane  intensity, 
the  aperture  plane  correlation  function  (Eq.  (24)),  and  we  have  also 
calculated  the  MTF  which  is  the  transform  of  the  PSF.  Next,  we  will 
use  the  MTF  and  the  degraded  image  spectrum  to  calculate  the  restored 
image  spectrum  (RIS)  IR(KX,K  )  using  the  approach  given  in  Eq.  (7). 


(29) 


[R(KX,Ky) 


Ii(KxiKy)  C(Kx»^y) 

H(Kx,Ky)  H(Kx,Ky) 


Substituting  Eqs. (15a), (15b) and  (15c)  into  Eqs.  (24)  and  (28)  and 
performing  the  division  indicated  in  Eq.  (29),  the  restored  image 
spectrum  is  calculated  in  terms  of  Kx  and  Ky,  the  spatial  frequency 
components;  a,,  b:  and  Cj,  the  object  phaseyfront  parameters;  and  Aalt 
Ab i  and  ac1?  the  reference  phase  perturbation  parameters.  The  notation 
used  for  this  RIS  given  be^ow  has  been  chosen  to  emphasize  the  eight 
independent  variables  discussed  above.  The  notation  to  be  used  is: 


(3°)  IR(Kx*Ky)  =  R(Kx«  ^y>ajjbi,C2,Aa|,Ab^,ACi)  , 


Actual  calculation  of  the  RIS  is  straightforward  but  tedious;  thus, 
it  is  carried  out  in  Appendix  B.  The  result  of  this  calculation  is 
that  the  RIS  has  the  mathematical  form: 


(31)  R(KX , K  ,aj jbj.CpAa^Abj jACj)  = 


1  + 


2Aa2Kx  +  AbjKy" 
2ai^x  +  bi^y_ 


Cos  {(2tt/K^  )^(2Km-|Kx| )  (2Aa1Kx+Ab1Ky)| 
-Sin|(2TT/Km2)J5(2Km-|Kx|)(2Aa1Kx+Ab1Ky)J.Cot|(27T/Km2)i5(2Km-|Kx|) 
•(2(a1+Aa1)Kx+(b1+Ab1)Ky)| 


2ac.Kv  +  Ab, Ky 

1  + - LI _ _ 
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2ciKy  +  biKx 


Co3-f(2n/Km?)^(2Kni-|Ky|)(2Ac1Ky  +  Ab^)}  -  Sin^ir/^) 
•ii(2Km-|Ky|)(2Ac1Ky  +  Ab1Kx)j-Cot|{2Tr/!<m2)5s(2Km-!Ky|) 
•(2(c1+AC1)Ky  +  (bj+AbjjKx)}]  . 


for  small  valuesof  Aa, ,  Ab  and  ac1  (i.e.,  much  less  than  1.0), 
Eq.  (31)  may  be  simplified  by  using  the  small  argument  approximations 
for  the  Sin  and  Cos  functions.  That  is,  for  small  values  of  the 
arguments,  we  have: 

Sin  a  -  a 

Cos  a  '  1 

where  a  may  be  the  argument  of  any  Sin  or  Cos  term  in  Eq.  (31).  Using 
these  approximations,  Eq.  (31)  is  greatly  simplified  and  may  be  re¬ 
written  as: 


(32 )  R (Kx ,Ky .ajjbj.CjjAaj, Ab  ^ ,Ac^ ) 


2Aa1Kx  +  AbjKy 

2aiK*  +  W 


•f1  -  {(2»/V>!'*(2VlKxl>(2A>1V''biKy)}Cot{(2'/l<m2)'s(2i'-™'|Kxl) 

2ac, Kv  +  Ab, K  " 

1  +  y  x 
2ciKy  +  bjKx 

L1  -{(2»/Km2)'s(2Kln-|KylH2«1V4blKx)}Cot{(2r/Km2)^(2K[n-|Ky|) 
.(2(c1+/sc1)Ky  +  (bj+abjlK,!]  . 


At  this  point  we  have  a  mathematical  relationship  which  describes 
the  RIS  in  terms  of  the  parameters  a  ,  b  ,  c,,  Aa,,  Ab,,  and  ac,  which 
are  associated  with  the  turbulence-irlduced  effects.  This  expression  is 
used  in  the  next  section  of  this  chapter  where  the  discrete  form  for  the 
RIS  is  derived;  it  is  also  used  in  Chapter  III  where  a  thorough  analysis 
of  the  RIS  is  conducted. 


•(2(a1+aa1)Kx  +  (bj+Ab, )Ky)}| 


17 


D*  The  Discrete  Restored  Image  Spectrum 

cno^<-The  Gxpressi01?  9  ven  in  Eq.  (32)  describes  the  restored  image 

pmnwS‘nna5-a^°5,t1nU'US  funct1oni  h°wever,  restoration  techniques 
2!  9.di9lt-  computers  require  only  a  finite  niznber  of  data 
samples,  taken  in  comp! nance  with  the  Nvquist  sampling  criterion,  to 
escribe  a  band-limited  function.  Since  we  are  interested  in  com- 
puter  processing  of  turbulence  degraded  images,  we  now  proceed  to  derive 
the  discrete  form  of  Eq.  (32)  which  we  assume  to  be  a  valid  repre¬ 
sentation  for  the  RIS.  M 

h  c,I?i -dnS^uetw  RIS,f°r.an  actual  imaging  system  would  be  obtained 
by  sampling  the  degraded  image  and  the  PSF  in  the  telescope  image 

The  rGsjlt;n9  two  sampled  data  arrays  are  then  read  into  a 
digital  computer  where  their  Fourier  transforms  are  calculated  using 
a  fast  Fourier  transform  (FFT)  algorithm.  (See  Appendix  C  for  a 
discussion  of  the  FFT.)  The  ratio  of  the  two  discrete  transforms 
U.e.,  image  spectrum  divided  by  MTF)  yields  the  discrete  RIS. 

The  sample  spacing  used  in  the  image  plane  and  the  quantitized 
spati3l  frequencies  are  next  calculated.  In  Eq.  (24)  the  spatial 
frequency  band  limit  for  the  image  intensity  was  found  to  be  2Km. 

For  this  study  we  assume  that  the  images  are  sampled  with  64  x  64 
sample  resolution.  Thus,  in  accordance  with  Nyquist's  criterion 
the  sample  spacing  must  be:  * 


(32a)  AX2  =  Ay2  Ti/2Km 

for  which  the  quantitized  spatial  frequency  is: 
(33b)  aKx  =  AKy  =  4y64  =  2tt/64ax2=  27r/64Ay 


(See  Appendix  C  for  details). 

.  is  assumed  for  the  model  used  here  that  a  telescope  input 

+?nncUi^  «tSally.11i;i£s  /-e  fre(iuency  content  of  the  intensity  func- 

nianS  *°  &i^Kx  and  63a*7-i  1,e*’  64aKx  =  64AKy  =  0);  however,  the  image 
plane  functions  are  sailed  with  64  x  64  resolution.  Thus,  the  image 

is  sampled  at  greater  than  the  Nyquist  rate,  which  eliminates  the  prob¬ 
lem  of  aliasing  in  the  spectrum. 

We  are  now  equipped  to  calculate  the  discrete  RIS  from  Eq.  (32) 
by  making  the  substitutions:  H  v  ‘ 


(34a) 

Kx  *  P4Kx 

(34b) 

Ky  - 
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where  p  and  q  -  0,  il»2»*’*31  and  AKX  and  AKV  are  the  discrete  spatial 
frequencies  given  in  Eq.  (33b).  y 

In  terms  of  the  "R"  function  notation  introduced  in  Eq.  (30)  the 
discrete  RIS  is  given  in  Eq.  (35)  below. 


(35) 


Discrete  RIS  =  R(p  aKv>  q  AlC, 

*  y 


a1*bi’c1»Aa1>Ab1>Ac1) 


The  restored  image  is  obtained  by  taking  the  inverse  FFT  of  the 
function  given  in  Eq.  (35).  It  is  the  quality  of  this  resulting  re¬ 
stored  image  that  we  discuss  in  Chapter  IV. 


It  is  worthwhile  to  show  the  relationship  between  the  discrete 
frequencies  in  the  RIS  passband  and  Km,  the  input  aperture  halfwidth, 
where  this  halfwidth  is  expressed  in  units  of  spatial  frequency  (see 
tq.  (lb)).  These  relationships  are  presented  in  Fig.  3  for  only  the 
.x  components.  The  Ky  components  are  exactly  the  same  since  a  square 
input  aperture  was  assumed.  M 


■RIS  PASS  BAND- 


zxxxxxxxxxxxxxxxx 


-INPUT  APERTURE  WIDTH- 


APERTURE 

-  HALF - 

WIDTH 


c|xxxxxxxxx  xx xxxxjcxxxxxxxxxxxxxxxjxxxxxxxxxxxxxxxxxl 
-3lAK-x«-2Kn,  — 15/z  Ax,  '-ATm  °  15  l/2  AKx=Km 

X«  DISCRETE  SPATIAL  FREQUENCY  VALUES 


31  A/Cx=2  Km 


Fig.  3.  Discrete  frequency  relationships. 
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The  RIS  passband  is  seen  to  be  twice  the  input  aperture  width  Km  which 
is  the  E-field  pass  band.  This  fact  results  from  convolving  the 
aperture  function  W(Kx,Ky)  with  itself  in  Eq.  (20). 

E.  Summary 

At  this  point  we  have  both  a  continuous  mathematical  model 
(Eq.  (32))  and  a  discrete  model  (Eq.  (35))  describing  the  restored 
image  spectrum  in  terms  of  the  parameters  ai,  bi,  ci,  Aai,  Abi,  and  aci 
which  describe  the  turbulence-induced  effects  in  the  received  E-field. 
The  goal  now  is  to  use  the  discrete  model,  in  conjunction  with  a  high* 
speed  digital  computer,  to  set  limits  on  the  ranges  over  which  Aai,  Abi, 
and  Ac,  may  vary  for  a  particular  choice  of  a,,  b,,  and  c,  such  that^ 
the  restored  image  (i.e.,  inverse  FFT  of  Eq.  (35))  is  a  "good"  approxi¬ 
mation^  the  actual  point  object.  The  details  of  this  "qoodness" 
criterion  are  discussed  in  Chapter  m,  section  E. 


CHAPTER  III 

ANALYSIS  OF  THE  RESTORED  IMAGE  SPECTRUM 


A.  Introducti 


on 


chapter  we  discuss  in  depth  the  concepts  associated  with 
thelr  qua1,t>'  and  the  restored  image  spectrum  (RIS) 

I™  «b,ch  fheV!’e  obtained.  First,  we  define  what  we  mean  by  an 
*  y-r^t0red,;m^  and  the  spectrum  of  this  image.  We  then  consider 
deaf  p  l  h“  wel1.the  « s.  9l*en  in  Eq.  (32)  of  Chlpter  II,  matches  The 

pofes  present  inm?herRK  d?h9ra?1n?,fa5‘?rs  tare  found  t0  ba  mathematical 
LrQ!fP^Sentin  i  RIS‘  Physical1y  this  theory  is  not  difficult  to 

RIS  PtIn1nrL^ht  9ross  departures  from  the  ideally  flat 

RIS.  In  order  to  facilitate  a  computer  analysis  of  these  Doles  a 

°f  e9uations  are  derived  relating  the  number  and  location  of 

(Eq  (14))S  t0  ^  parameters  Ascribing  the  E-field  phase  expansion 

described^rchan^r'ri?6  effectl' veness  of  the  restoration  procedure 
-  5  Chapter  1 »  ]  t  is  necessary  to  describe  a  criterion  with 

chaDter°wi th^3  ^stored  ^age  quality.  Therefore,  we  conclude  this 

restored  l^age  q^’"  °f  three  paramet-rs  chosen  *«ribe  the 


B. 


Ds  i  i  n  i  ti  on  of  the  Ideally  Restored  Image 


given  for  the 
scussion  of  the 


In  this  section  the  ideally  restored  image  is  defined  for  the 
degraded  aperture  plane  E-field  model  presented  in  Chapter  I  The 
mathematical  expression  for  the  ideally  restored  image  is  -<« 
continuous  and  discrete  cases.  And  finally,  a  brief  disc 
non-ideally  restored  image  is  presented. 

The  ideally  restored  image  is  defined  in  terms  of  the  equations 

DlSnTE1f?olHeH^trt’“ie?“  E-field  wvefronts.  lhese  aperture 

p  ane  E-field  distributions  for  the  object  point  and  reference  Doint 

are  reproduced  below  in  Eqs .  (36a)  and  (36b),  respective!™  P 


(36a)  Ea(Kx,Ky)  =  e 


j(dj/k)2(aKx2+bKxKy+cKy2) 


(36b)  Ea(Kx,Ky)  = 


j(di/k)2((a+Aa)Kx2+(b+Ab)KxKy+(c+Ac)Kv2) 


^hlhLSPeCia1^aS!uin  whlch  Aa  "  Ab  =  Ac  =  0,  both  Eqs.  (36a)  and 
(36b)  have  exactly  the  same  mathematical  form.  This  means  that  both 

?hrQunHe^P^rtj;1edt1'n+exaCut1y  the  saTO  manner  **  the^  propagated 
la  lent  atmosphere.  Thus,  using  the  MTF  obtained  from 

Eq.  (36b)  to  restore  the  degraded  image  characterized  by  Eq.  (36a),  we 
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define  the  ideally  restored  image  as  the  image  obtained  when  Aa  =  Ab  = 
Ac  =  0.  This  definition,  an  extension  of  Harris's  [9]  work,  does  not 
restrict  the  values  of  a,  b,  and  c.  Physically,  this  means  that  as 
long  as  both  the  reference  wave  and  the  object  wave  experience  the  same 
turbulence  effects,  as  characterized  by  a,  b,  and  c,  the  best  possible 
restored  image  is  obtained. 

Applying  the  condition  that  Aa  =  Ab  =  Ac  =  0  (equivalently 
Aax  =  Ab  =  Acx  =  0)  to  the  continuous  RIS  in  Eq.  (32)  and  taking  the 
inverse  Fourier  transform  of  this  function,  an  expression  is  obtained 
for  the  ideally  restored  image. 


1  ((  j (K„x  +K  y  ) 

(37a)  IR(Xl,y.)  =  -i— JJ  R(Kx,Ky ^ ,bj  ,Cl ,0,0,0)e  y  dK 


(37b) 


(37c) 


1 

^n; 

rr 

(2tt  )2 

jj 

"^m 

1 

rr 

<2„)2 

JJ 

®m\2 

Sin 

xdKy 


1  (( 

-prjj  le  y  dKxdK} 


2«mx1  2^! 


Note  that  the  result  obtained  in  Eq.  (37c)  is  not  the  ideal  point 
object,  5(xi,yi);  however,  in  the  limit  as  Kjn  ®  ,  Eq.(37c)  dees 
approach  this  ideal  result,  a  fact  well  known  from  diffraction  theory. 
Thus,  the  system  band  limit  2Km  has  imposed  an  upper  bound  on  the  re¬ 
stored  image  quality  even  in  the  ideal  case. 

When  the  discrete  case  (Eq.  (35))  is  considered  for  Aa  =  Ab  =  Ac 
=  0  (equivalently  Aa:  =  Abx  =  Acx  =  0),  we  must  recall  that  the  FFT 
algorithm  represents  the  function  being  transformed  as  a  discrete 
Fourier  series  (see  Appendix  C) .  The  results  indicate  that  when 
R(p  AKX,  q  AK  ,  aj,  b^  cx,  0,  0,  0),  which  equals  1.0  at  all  discrete 
frequencies,  Is  inverse  Fourier  transformed;  the  resulting  function  is 
the  discrete  analogy  to  the  Dirac  delta  function.  For  this  case,  the 
discrete  restored  image  is: 


(38)  IR(p  AXj ,q  AXj)  =  6 (p  AXj.q  Ay^ 
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c. 


where  p  and  q  are  integers  that  specify  the  location  of  the  restored 
image  in  the  6  x  64  sample  restored  image  array. 

Neither  Eq.  (37c)  nor  (38)  is  the  ideal  ^xi.y,)  point  obie-t  but 
ftfthe  band°lin,??p7?retent  the  beSt  poss1b,e  aPP^™ation  of  f(x.,yi) 

.WMS  STS  ^  «"  ^ 

2^n??S1iii?£?hj°tnoer 

need  fo?  ^  i  1*“^  ’“S'-  ™s  fact  ™Phasi2es  the 
neeo  tor  a  analysis  of  the  degrading  effects  of  finite  Aa,,  Ab:  and 

Analysis  of  the  Restored  Image  Spectrum  (RIS) 
lentl^Ar^T  Llf suggeJ.^d  that  when  4a,  Ab  and  ac  (or  equiva- 

swS?:Sat«|®^«S: 

5.*r»s;‘r,:K‘&„,„,K'B  rib?£B 

The  restored  image  obtaineo  when  Aa,,  Abi  and  Ac,  are  finite.  ic 

Wh™nLT=tedb*0='>nc"^1deR1(KbyKreaU™1'nc9  ‘n  o‘hmR!S  9iV6'!  in  Eq'  (3?)' 

?rLfft^ndCiy?eldSthiKetdeSifW’PrSb"b:  W^.’Sta’inwn. 

H  u*  yie1ds  Jhe  1  dea^y  restored  image  calculated  in  Eq  (37 c) 

™Tv^g’thtencot4aiunc  l„„a"dn\Cl-are,f,'n^-  the  RIS  i"clu^  teL  in- 
conditloJ  ocSrs  (?  e  a  pX  eldstsl  L'  Cf°‘  '  (n  1nt*9er) 

greatly  exceed  the  rieciJLPw!i  f  the  frequency  component  will 

“eai'r^to^LV™5'0™  °f  ?his  *3ra*dreRISt°„l1V?a??  the  & 

.JI,?.l“nt,n“,'5  f0™  for  the  RIS  is  “sed  this  section-  however 

slisSSf™  ssu;.!i,?i  -*• 

the  function  Z(Kx,Ky),  assumes  the  value  Z  =  nn  (n  =  integer). 
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(39)  R ( Kx , Ky , a j  .bj.Cj ,Aaj  ,Ab1,Ac1 ) 


ZAa^  +AbjKy 
2ajKx  +  bjKy 


1  - 


(2ir/Kffl?)!S(2!(m-|Kx|)(2Aa1Ky+Ab1Kv)|-Cotj(2^/Km)12(2Km-  K^) 


{- 


•(2(a1+Aa1)Kx  +  (bj+Abj  )Ky)[ 


1  + 


2a  c  j  Ky+  Ab  j  Kx 


2c1Ky+b1Kx 


1  -{(27T/Km2)1^(2,<in-|Ky  |  )(2AC1Ky+Ab1Kx)|cot|(2TT/Km2)Js(2Km->Ky  | ) 
•(2(c1+AC1)Ky  +  (bj+Abj )KX)||  . 


Since  both  square-bracketed  terms  involving  the  Cot  function  have 
the  same  mathematical  form,  we  need  only  analyze  one  such  term.  The 
other  term  will  follow  with  the  proper  substitution  of  wavefront  para¬ 
meters.  Thus,  we  need  only  analyze  the  term  involving  the  a  and  b 
parameters.  The  Cot  argument,  Z(KX,KV),  for  this  term  is  given  in 
Eqs.  (40a)  and  (40b)  below.  y 


(40a)  Z(Kx,Ky)  =  (2tt7 Km2 )^(2Km- [ Kx | ) ( 2apKx  +  bpKy) 

(40b)  -  (’T/Km2)^-2apKx|Kx|-bp|Kx|Ky  +  4apKmKx  +  2bpKmKyJ 

The  parameters  ap,  bp  and  cp  are  defined  below  as: 

(41a)  ap  =  aj  +  Aaj 

(41b)  bp  =  bj  +  Ab  j 

(41c)  cp  =  c1  +  ACj 

Eq.  (40b)  is  an  odd  function  of  Kx  and  Ky;  that  is, 

Z(Kx,Ky)  =  -Z(-Kx,-Ky). 
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Thus,  the  analysis  may  be  further  limited  to  the  half  plane  defined  by: 

Kx>0  ,  ■2Km<  Ky<  2Km  . 

The  analysis  for  the  second  half  plane  is  obtained  by  substituting 
-Kx  for  Kx  and  -Ky  for  Ky  into  the  equations  derived  below  and  then 
multiplying  by  -1.0. 

As  stated  previously,  the  loci  of  poles  are  determined  by  the 
condition  Z(Kx,Ky)  =  nTr  (n  integer).  The  maximum  number  of  such  pole 
loci  is  obviously  proportional  to  the  maximum  range  of  Z(Kx,Ky). 
Knowledge  of  this  range,  in  terms  of  ap  and  bp,  affords  control  over  the 
number  of  pole  loci.  This  facili tates  a  computer  analysis  of  the 
effects  of  these  parameters  on  the  restored  image  quality. 

In  order  to  determine  this  maximum  range,  we  must  first  generate 
an  equation  in  terms  of  ap  and  Up  and  determine  the  maximum  and  minimum 
values  assumed  by  Z(Kx,Ky;.  We  will  find  that  this  analysis  is  divided 
into  two  cases  depending  on  the  relative  magnitudes  of  ap  and  bp. 

Using  the  first  and  second  derivatives  of  Z  with  respect  to  Ky  to 
determine  maxima  or  minima  with  respect  to  Ky,  we  find  that: 


(42a) 

and 


SlL—  =  o 

3Ky 


Kx  =  2K,, 


(42b) 


32Z  = 
3Ky2 


Substituting  Eq.  (42a)  into  Eq.  (40b)  we  find  that  Z(Kx,Ky)  =  0  for  all 
Ky.  Since  Eq.  (42b),  the  second  derivative  test  for  maxima  or  minima 
fails,  we  deduce  that  there  are  no  maxima  or  minima  along  the  line 
Kx  =  2  Km . 

From  the  first  and  second  derivatives  of  Z  with  respect  to  Kx  we 
find,  with  respect  to  Kx,  that: 


(43a) 


aZ_  _ 
3KV 


4aP 

Ky  =  TT  (KnfKx) 


and 


(43b) 


32Z 

3KX2 


4a. 
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a/h’  !hn  5eCOnd  ^erivative  test  indicates  that  maxima  do  exist  for 
TfnsPequation  a,°"9  the  ,1ne  Ky=(4ap/bp) (Km-Kx) 


K1  •’bf'  <  *m-  *,  > 


-(TT) **■ 


-2k, 


LINE  OF  MAXIMA  FOR-P>( 

bP 

LINE  OF  MINIMA  FOR  ^  <  C 


Fig.  4.  Maxima  and  minima  of  Z(Kx,Ky). 

Evaluating  Eq.  (40b)  along  the  line  Ky  =  (4ap/bp)(VKx),  we  have: 

(44)  Z<Ky’KJ'  *  "bf  (K"'-Kx))  =  (*«n,2)(2apKx2  -  8apK„,Kx  ♦  8apKm")  , 

2^  onA  &»*•«>. 

from  ""  Z  Ve™  Variatio"  w1th  Kx  -  0,  we  find 


Fig.  5.  l  evaluated  along  line  of  maxima  or  minima  vs.  K 
(45)  Z(0,Ky)  =  („/Kmt)  ZbpK^Ky  ,  1 

which  has  a  maximum  at  K  =  2K„,  as  seen  in  Fig.  6.  Thus,  when  b„  >  2an, 

hS  ,n  *Tv.Va!-e  a‘  !°-2Km).  Since  Z  has  been  showS  to  P 

De  an  odd  function  of  Ky,  the  minimum  must  occur  at  the  point  (0,-2Km) 

as  shown  in  Fig  6.  Thtis,  the  maximum  range  of  Z,  the  quantity  we  are 
interested  in,  for  bp  >  2ap  is: 


^max  “  ^min  =  “  ( ~^t>pTT )  =  8bpiT 

The  parameter  bg  in  Eq.  (46)  determines  the  maximum  number  of  pole  loc1’ 
contributed  to  the  restored  image  spectrum.  P 

from  Eo^mRxI  fJa’c?  situation  exists.  The  results  obtained 

om  Eq .  (43a )  ana  Fig.  4  indicate  that  the  maximum  value  of  Z  with 
respect  to  Kx  occurs  at  the  intersection  of  the  lines  Ky  =  2Km  and 


27 


Fig.  6.  Zmax  and  Zml-n  as  a  function  of  Ky. 


K„  =  (4ap/bp)(Krn-Kx)  or,  by  substituting  Ky 


y  '  P'  p/' 'Tn  ''x/  u 
(Km-Kx),  at  the  point: 


2Km  into  Ky  =  (4ap/bp) 


(47a)  Kx  =  Km(l  -  -£-) 

cap 

(47b)  Ky-ZK,  . 


This  result,  which  is  valid  for  -2ap<  bp<  2ab,  is  shown  in  Fig.  7. 

Consider  a  (Ky,Z)  plane  cutting  through  Fig.  7  along  the  line 
Kx  =  Kxg .  By  substituting  Kx  =  Kx  into  Eq.  (40b)  and  regrouping  terms, 
we  obtain:  u 

(48)  Z  (KXq ,Ky )  =  (tt/Kjd2  )  |^(2bpKm-bpKXo)Ky+4apKniKXo-2apKx2  , 
which  is  graphed  in  Fig.  8. 


2  SLOPE  s  (2bp  Ktn-  tnKt0 )  */<* 


t  °p*m  *io“20pK,0l  /K 


I  t 

2op  **o_^°p*m*io 


2  min*  /k  ^  [j~  (2bp  Km  -  bp  Kt0 )  2 
1  +4ap^m^o“2oP/C/o] 


Fig.  8.  Zmax  and  Zm-jn  as  a  function  of  Ky. 

In  Fig.  8  it  is  seen  that  the  maximum  value  of  Z,  with  respect  to 
Ky,  occurs  along  the  line  K  =  2Km.  In  Fig.  7,  it  was  demonstrated 
tnat,  for  -2ap<  bp<  '2ap,  the  maximum  value  of  Z  with  respect  to  Kv 
occurs  along  ?he  iTne  K/=  Km(l  -  (bR/2ap)).  Evaluating  Zmax  (from* 
Fig.  8)  at  the  point  (Kx  =  Km( ] -(bp/2ap ) ? , Ky  =  2Km) ,  we  fiml? 

T  b  2" 

^max  =  2ap  +  2bp  +  -2—  tt 

2ap_ 

Since  the  Z  axis  intercept  in  Fig.  8  will  always  be  greater  than 
or  equal  to  zero  for  0  <  Kx  <  2K|j,,  the  minimum  value  of  Z  will  always 
occur  (see  Fig.  6)  at  the  point  fkx  =  0,  Ky  =  -2Km),  thus: 


(50)  Zmin  =  -4tyr 

Using  Eqs.  (49)  and  (50),  the  maximum  range  for  Z,  when  -2a„<  b  <  2a 
is  calculated  to  be:  P  P  P’ 


^max  "  ^min  =  (^ap  +  6bp  +  ^  )  n 

P 

In  this  case  the  magnitudes  of  both  the  parameters  aD  and  bD  contribute 
to  the  total  number  of  pole  loci  present  in  the  RIS.  Dividing  the 
maximum  ranges  for  Z  found  in  Eqs.  (46)  and  (51)  by  tt,  and  then  sub¬ 
stituting 


we  find  that 


ap  = 

bp  =  bj  +  Abj  , 

the  maximum  number  of  pole  loci  Mp  (M 


integer) 


is: 


(52a)  for  2(aJ  +  Aa^  <  bj  +  Abj  , 
(f'2b)  Mp  =  8lbi  +  AbJ 


and  for 


(53a)  -2(a1  +  Aaj)<  bJ  +  Abj  <  2(aT  +  Aaj)  , 

I  (b  +  Ab  )2 

(53b)  Mp  =  1 2 (a  +  AaT )  +  6 (b  +  Ab  )  +  — l - 3: 

1  *  A3,) 

The  gross  characteristics  of  the  pole  loci  in  the  RIS  are  determined  in 
Eqs. (52b)  and  (53b)  by  the  parameters  a,  and  b,,  which  are  a  measure  of 
the  wavefront  distortion  at  the  edge  of  the  input  aperture  in  fractions 
of  a  wavelength;  however,  the  perturbation  parameters  Aa,  and  Ab,  con¬ 
tribute  in  two  important  ways.  First,  they  may  increase  or  decrease 
the  value  of  Mp  by  a  sufficient  magnitude  to  actually  introduce  or 
remove  a  loci  of  poles.  Second,  regardless  of  whether  they  do  introduce 
or  remove  a  loci  of  poles,  the  perturbation  parameters  will  always  tend 
to  shift  the  location  of  the  existing  loci  of  the  poles  in  the  KXKV 
plane.  This  second  point  is  important  in  terms  of  the  discrete  caBe  in 
which  a  small  change  in  Aaj  or  Abj  may  result  in  a  discrete  frequency 
component  being  evaluated  precisely  at  a  pole.  The  implications  of 
this  event  are  covered  in  Chapter  IV,  section  B. 
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analv?i!  of^hi^rJ  sho“ld  be  covered  before  completing  this 

analysis  of  the  RIS.  This  point  involves  the  term 


(54)  aZ  =  7r/Km2[-2Aa1Kx|Kx|-Ab1|Kx|Ky+4Aa1KmKx+2Ab1KmKy] 

in'fLni^r!165  ti-  C°r  fu"ction  in  Ec!-  (39)«  It  was  stated  previously 
thaPl^n  Th*  Sec^1on  C;  thet  the  ranges  for  Aa,  and  Abj  are  much  less 
tional’Sn  slJce.Afand  (^Z/dKx)  and  (9AZ/%)  are  all  propor- 

Z  is  siali^nH  iUrb*  •  Param?ters>  this  implies  that  the  function 
Z  is  small  and  slowly  varying  in  both  the  Kx  and  Kv  directions.  This 

iCatlH6  th6n  tends  t0  rnodu1ate  the  Cot  function  in  such  a 

Ln  Ke  Hce  th!  undr  the  Cot  faction.  In  the  next  chapter 
this  will  be  demonstrated  to  be  particularly  important  near  RIS  poles. 

The  most  significant  results  obtained  from  the  RIS  analysis  are 

thevl1rel'at2Sthi52b  -and  *  These  epuations  are  emphasized  because 

oStheslzeS  tn  h/f lmUm  nU?6r  of  RIS  P0le  1oci’  whictl  have  been  hy- 
.  ’z  d  to  have  a  significant  bearing  on  the  restored  image  quality, 

(see  Eos  ?3fiaHnHS?tnideS^1be  th?  turbulence  degraded  input  waves 
In  6  l-  d  P6bP*  These  relationships  will  be  used  in 

Chapter  IV,  section  B,  where  actual  image  restoration  is  performed 

on  three  classes  of  degraded  images. 

D.  The  Defocused  Image  Case 

fnr  P°j-’  a  particular  class  of  degraded  images  is  selected 

for  analysis.  Tms  class  of  images  is  distinguished  by  the  absence  of 

the  coupling  between  Kx  and  Ky  in  Eqs.  (36a)  and  (36b)!  That  is  wp 
set  the  parameters  b  -  Ab  =  0  (or  equi valently  bj  =  Ab,  =  0)  in  all 

tionriPrnH^-V1H9-th^PaVefr0nt  parareters.  A  review  of  these  equa¬ 
tions  is  contained  in  this  section. 

rerordinn^nf^ Hof  this  case  (i.e.,  bj  =  Abj  =  0)  corresponds  to  the 
is  vlr!  w  1-ageS  !n  the  telescope  image  plane.  This  case 

the  RTS  an hPJ  wher\?Ie  ’S  studying  the  relationship  between  poles  in 
the  RIS  and  the  quality  of  the  resulting  restored  image.  It  is  useful 
because  the  associated  mathematics  are  simpler  than  those  for  the  case 

theVPand  thpPP  t^H  ^  PAbl  bUt  yet  l0ci  of  poles  are  present  in 
the  RIS  and  they  do  tend  to  degrade  the  restored  image  quality. 

rhantPPu  PI1??!"9  dlscussion,  the  general  equations  derived  in 

S?mEv1  leSfv  1  and  ?u\atei  cross-coupling  terms  b  and  Ab 

\or  equivalently  bj  and  Abj)  set  equal  to  zero. 

The  aperture  plane  E-field  distribution  taken  from  Eq  (36a)  is 
rewritten  for  this  special  case  in  Eq.  (55)  below. 
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(55) 


Ea(Kx,Ky ) 


j  (di/k)°(aKx°  +  bK  2) 
=  e  J 


Assuming  this  phase  distribution,  the  degraded  image  spectrum  (aperture 
plane  correlation  function,  C(Kx,Ky),  and  the  MTF,  H(KX,KV),  have  the 
following  forms  (the  general  equations  are  (24)  and  (28) )y 


(56)  C(KX,K  )  = 


(2Km-|Kx|) 


Si  n  ( d-j  /k ) 2  (2  Km- 1  Kx  | )  aKx 
(d-j/k)2(2Km-  |Kx|)aKx 


(2Km-|Kj) 


(d-j/k )2(2Km- 1 Ky  |  )cKy 


(57)  H(Kx,Ky) 


|(2Km-|Kx|) 


S^n(di/k)2(2Km-|Kx|)(a  +  Aa)Kx 
(d-j/ k ) 2 ( 2Km-  |KX  [ )  (a  +  Aa)Kx 


(2Km-|Kvl>: 


Sin(di/k)2(2Km-|Kv|)(c  +Ac)K 


( di /k  N2 (2  Km- I Ky | ) ( c  +  Ac)Ky 


Note  that  Eqs.  (56)  and  (57)  are  separable  into  two  factors,  one  in¬ 
volving  only  Kx  and  the  other  involving  only  Kv.  Thus  the  behavior 
of  each  factor  is  graphed  in  a  single  figure  (fig.  9),  where  K  repre- 
sents  either  Kx  or  Ky  and<Y  represents  a  or  c.  The  total  function  is 
then  obtained  by  multiplying  two  such  diagrams  together  at  each  point 
in  the  KxKy  plane. 

(frorn^q  "(32 J)^  Spectrum  (RIS>  for  the  special  case  becomes 

(58)  R(Kx,Ky,a.,0,c.,4a1,0,ac1)  =  C_^y±_  = 

H(Kx,Ky) 


<2Knf  ! Kx I  HKxjcot {(2-r/V )(2K„- 1 Kx|  )apKx|) 
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Mg.  9.  Graph  of  (2Km-|K|)  Sinc[(di/k)2  (2Km- |K|  )YK], 

(1+^p)  (1'{(2’,/K|":’)(2Km'|Kyl,4cil<y}Cot{<2"/l<ni2)(2Kn|-|Kyl)cpKy| 

Note  that  each  factor  in  Eq.  (58)  involves  only  Kx  or  Kv;  thus,  the 
overall  inverse  transform  is  the  product  of  two  one-din&isional  inverse 
transforms. 

a  /rI!ik\dl-Sfrere  1S  obtained  by  substituting  Eqs.(34a) 

and  (34b)  into  Eq.  (58)  to  yield:  a  H  v  ' 

(59)  R(p  aKx,  q  AKy,  a^,  0,  Cj,  Aa^,  0,  Ac^) 

rnc  J5liUnibe//?Iuf1eJoci  for  this  special  case  is  obtained  from 
Eqs.  (52b)  and  (53b).  Obviously  bp  =  0  satisfies  Eq.  (53a);  thus  the 


number  of  pole  loci  is 

(6°a)  Mpx  =  2{3l  +  Aa j )  =  2 |ap |  . 

nf‘ni?0a^indiCutuSi^ha1t  the/x  term  in  Eq.  (58)  contributes  2an  loci 
of  poles  in  each  half  plane  (i.e.,  Kx  <  0  and  Kx  >  0).  Simi 1  a rP results 

rVlV  '  Kj,nle™  j"  Eq-  <5S)-  ^  two  K  half  planes 

1,e*’  y  —  9  an^  :L  0)  each  contribute  y 

(60b)  Mpy  =  2(ci  +  aCl)  =  2 |cp | 

loci  of  poles.  The  number  of  pole  loci  equals  minimum  integer  value 
or  tqs.  loUaj  and  (60b). 

solving theCteqoa{?oCn:°f  P°'eS  ^  den,onstrated  *°  ba  straight  lines  by 

(61)  Z(Ky,Ky)  =  nn 

tQ?n  e&Vn^tMc  C°?  TST"1  fr0m  Eq*  (58)  and  n  Is  an  integer. 
by  0  <  k  ?  Ik  aiSd  ""l1*!*  the  Kx  half  plane  defined 

ditions  we  obtain  1,1  'Ky<‘  2Km *  So  V1ng  E0-  (61)  for  these  con- 

(62a)  Z(KX ,Ky )  =  -  K  2  +  K  s  nr 

K  2  K  X 


which  reduces  to 


2ap-n 


(62b)  Kx  =  K,,,  ±  Kn 


Eq.  (62b)  is  indeed  the  equation  for  two  straight  lines  located  eoual 
distances  on  either  side  of  the  line  K  =K  tJ  l  L  h  q  • 

Kx  half  plane  and  the  two  K„  half  p?aneXs  Tillow  f  ™  s  milar  cad-9 
tions.  An  overall  view  of  these  pole  loci  is  presented  n  Fia  in  Th* 

ariserSwhenxthedCotVaSUS  ^  faphs  ?,ear,y  111“strate  that  the'po?e  loci 
arise  wnen  the  Cot  argument  Z  =  m  (n  =  1  for  Fig.  10).  K 

valued •  t5£  w-  aB  is1  Ceremented  by  a  small  amount  to  a  new  larger 
value  ap  ,  the  loci  of  poles  shift  in  the  RIS.  9  r 
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RESTOREO  IMAGE  SPECTRl 
Kt  K,  PLANE 


2*m 

i  n  n — 

I  P“‘ " 

At  this  point,  all  the  significant  equations  from  Chapter  II  and 
Chapter  IIIA,  B  and  C  have  been  rewritten  for  the  special  case  in  which 
there  is  no  cross  coupling  between  Kx  and  Ky  (i.e.,  bj  =  Ab i  =  0). 

Note  that  Fig.  10,  which  illustrates  the  location  and  shifting  of  pole 
loci,  is  a  very  important  graph  and  it  is  referred  to  several  times  in 
Chapter  IV  where  actual  image  restoration  of  degraded  images  having 
zero  cross  coupling  between  K  and  K  (i.e.,  defocused  image  case)  is 
discussed.  x  J 


E-  Restored  Image  Quality  Criterion 

tho  s®ctl0!J.ls  devoted  to  a  discussion  of  the  criterion  by  which 

thpm^it0re?  q+Lbty  1S  Judged.  Among  the  authors  who  have  addressed 
themselves  to  this  problem  have  been  Linfoot  f 16 |  and  Barry  [17|.  They 

a ve  emphasized  that  it  is  desirable  for  such  a  criterion  to  reflect  a 

auali F6  r^er  than  subjective  measure  of  the  restored  image 

ty\  For  thls  reason,  three  parameters,  the  restored  Strehl  ratio 

r!St0red  )mge,  integral  scale  M  and  restored  image 
“  qua?;rror  {cms )» were  selected  to  describe  the  restored  ?mage 
quality.  A  discussion  of  each  of  these  parameters  and  the’>  optimum 
values  comprises  the  material  discussed  in  this  section 


(63) 


The  Strehl  ratio  is  defined  as: 

Ii (0,0) 
s  =  - 

ll(0,0) 


where  ^(0,0)  is  the  point  spread  function  (PSF)  of  the  turbulence 
degraded  imaging  system  and  I j(0,0)  is  the  PSF  for  the  ideal  turbulence 
free  system.  However,  for  this  study  in  which  we  are  interested  in 
restored  image  quality,  we  define  a  restored  Strehl  ratio  Is: 


(64) 


Sr 


IR(0,0) 

Ii  (0.0) 


iS  t!!e  ima9e  Obtained  from  the  RIS  and  It(0,0) 

is  the  ideally  restored  image.  For  the  optimum  case  (with  a  point 
source  input)  the  ideally  restored  image  has  the  value  It (0,0)  =  1  0 
S  S  condition  irto  Eq.  (64),  we  see  that  tie  restored 
ra*'°  Sriuls  rm;rely  thc  restored  image  value  at  the  coordinate 
Sr =  Tr(oIo) "=  iTSUS>  f°r  the  ideally  ^stored  image  (see  Eq.  (38)), 

...  The  res  to  re  d  Strehl  ratio  yields  both  spatial  domain  (x1>yi)  and 
a^  f^uency(|<x»,<y). information  about  the  restored  image!  In 
terms  of  the  spatial  domain,  Sr  is  the  ratio  of  the  central  intensity 
peaks  of  the  restored  image  and  the  ideal  image.  For  the  spatial 
frequency  interpretation,  we  recall  that,  In  terms  of  Fourier 
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transforms,  IR(0,0)  and  I j (0,0)  may  be  written  as: 


(65a) 


and 


IR(0,0) 


(65b) 


Ir(0,0) 


1.0  dKx  dKy 


Eqs .  (65a)  and  (65o)  are  the  average  value  of  their  respective  spatial 
frequency  spectrum.  This  means  that  the  restored  Strehl  ratio  may  also 
be  thought  of  as  the  RIS  average  value  since  the  average  value  of  the 
ideally  restored  image  spectrum  is  1.0  (all  frequency  components  in  the 
spectrum  of  the  ideally  restored  image  are  equal  to  1.0).  Thus,  the 
restored  Strehl  ratio,  which  has  the  optimum  value  of  1.0,  gives  us  a 
measure  of  how  closely  the  restored  image  and  the  RIS  approach  their 
ideal  values. 


The  integral  scale  r0  is  defined  as: 


(66) 


r0  =  + 


63  63 

l  l  Ir(P  Ax,  q  Ay) 
Pgl  q=i _ 


TT  IR(0,0) 


where  r0  is  defined  such  that  the  volume  of  a  cylinder  with  height 
lp(0,0)  (restored  image  value  at  the  origin)  and  radius  rn  equals  the 
volume  under  the  restored  image  intensity  function  IpU^y,).  This 
definition  is  illustrated  in  Fig.  11.  The  parameter  r0  gives  a  measure 
of  the  spread  of  the  restored  image.  The  optimum  value  of  r0  is 
r0  =  1/ At  which  occurs  when  the  integral  in  Eq.  (66)  equals  IR(0,0). 

If  Ecl-  (65)  is  divided  by  1/ /if,  a  normalized  integral  scale  is 
defined  as 


(67) 


J 


for  which  the  optimum  value  of  r  is  r  =  1.0.  It  is  the  normalized  in¬ 
tegral  scale  that  will  be  used  in  this  study  to  express  the  spread  of 
the  restored  image. 
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Fig.  11.  Diagram  illustrating  the  definition  of  r 
the  integral  scale. 


The  final  restored  image  quality  parameter  is  the  restored  i 
mean  squared  error  (ems).  This  parameter  is  defined  as: 

,  63  63 

<«>  'ms  -  <5r)2  l  l 

*  P=0  q=0 

where  I R  and  I j  are  the  restored  and  ideal  images,  respectively. 

The  parameter  ems  gives  a  measure  of  the  mean  squared  error 
between  the  actual  restored  image  and  the  ideally  restored  image. 
Obviously  the  optimum  value  for  is  zero. 


|lR(p  ax,  q  Ay)  -  I  j  (p  ax,  q  Ay)j2 


Next  in  the  discussion  of  the  restoration  parameters  is  the  appli¬ 
cation  of  the  computer  processing  techniques  to  the  ideally  restored 
images  defined  by  Aa,  =  Ac,  =  0  and  finite  a,  and  c,  (see  Chapter  III-B). 
The  resulting  restoration  parameters  calculated  using  Eqs.  (64),  (67) 
and  (68)  are  defined  as  bench  marks  against  which  the  parameters  calcu¬ 
lated  for  images  having  finite  valued  Aaj  and  Aq  may  be  compared. 

After  numerous  computer  runs  using  various  values  of  a,  and  c, 

(Aa  =  Ac,  =  0),  i t  was  verified  that  the  bench-mark  values  are  inde¬ 
pendent  of  the  values  selected  for  a,  and  Cl . 

Table  1  lists  both  the  theoretical  and  the  actual  computer 
calculated  values  for  each  parameter. 

TABLE  1 


Parameter 

Calculated  Value 

Theoretical  Value 

(Aai  =  ac,  =  0) 

Sr 

1.000 

1.000 

r 

1.000 

1.000 

cms 

0.000 

8.553  x  10'10 

The  discrepancies  between  the  theoretical  and  calculated 
values  of  ems  represent  inherent  computer  processing  errors;  however, 
these  errors  are  sufficiently  small  to  be  insignificant.  Thus,  the 
calculated  values  are  accepted  as  the  bench-mark  parameters. 

The  rang  j  for  Sr,  r,  and  which  define  a  successfully  restored 
image  must  be  specified.  A  successfully  restored  image  will  be  defined 
in  this  study  as  an  image  with  the  following  properties; 


(69a)  0.5  <  5r  <2.0 
(69b)  0.5 <  r  <  2.0 
(69c)  <  0.1 


The  range  for  Sr  allows  a  ±3  db  deviation  in  the  restored  image  central 
peak  intensity.  The  spread  of  the  restored  image,  as  expressed  by  r 
may  vary  by  a  ratio  of  2:1  while  the  mean  squared  error  must  be  less’ 
than  10  >.  These  values  for  the  quality  parameters  were  selected  on  the 
basis  of  preliminary  computer  processed  images.  Obviously  either  more 
stringent  or  more  relaxed  requirements  could  have  been  imposed;  how- 
ever.  this  particular  choice  serves  as  a  reasonable  starting  point  for 
defining  a  successfully  restored  image. 


40 


This  completes  the  discussion  of  the  restored  image  quality 
criterion.  The  quality  parameters,  restored  Strehl  ratio  (Sr),  integral 
scale  (r),  and  mean  squared  error  (c^),  were  defined  and  the  limits 
within  which  successful  image  restoration  occurs  were  specified  for 
these  parameters. 


This  chapter  has  been  primarily  devoted  to  the  discussion  of 
restored  image  spectrum  ( RIS)  and  the  restriction  of  this  general 
analysis  for  the  defocused  case  in  which  there  is  no  cross  coupling 
between  the  Kx  and  Ky  variables.  The  most  significant  result  of  this 
analysis  is  a  series^of  equations  (Eqs .  (52b),  (53b),  (60a)  and  (60b)) 
relating  the  number  of  pole  loci  in  the  RIS  to  the  wavefront  Darameters 

v  %  and  y 

A  discussion  of  the  ideally  restored  image  and  the  restored  image 
quality  criterion,  which  determines  whether  a  processed  image  may  be 
classed  as  a  successfully  restored  image,  were  also  discussed  in  this 
chapter. 


CHAPTER  IV 
RESULTS 


A.  Introduction 

In  this  chapter,  three  classes  of  turbulence  degraded  images  are 
restored  using  a  computer  processing  technique  which  involves  the 
Fourier  transform  method  of  image  restoration.  (This  method  was  dis¬ 
cussed  in  Chapter  I.)  Each  class  of  degraded  images  is  defined  by  the 
number  and  location  of  the  pole  loci  present  in  the  RIS.  The  success 
or  failure  of  the  attempted  restorations  is  then  discussed  in  terms  of 
the  restored  image  quality  criterion  (see  Chapter  III,  section  E).  Graphs 
which  present  the  behavior  of  the  restored  image  quality  parameter  as  a 
function  of  the  wavefront  perturbation  parameters  Aa  1  and  aCj  are 
included.  The  allowable  ranges  for  Aaj  and  ACj  over  which  successful 
image  restoration  does  occur  are  defined  with  reference  to  these  graphs. 
These  limits  on  Aa!  and  acj  are  the  primary  goal  of  this  study  and 
they  will  be  used  to  define  the  ranges  of  successful  image  restoration. 

B.  Restoration  Results 

It  has  been  suggested  in  Chapter  III,  section  C,  that  the  presence 
of  poles  in  the  restored  image  spectrum  (RIS),  as  determined  by  the 
parameters  ap(=a!  +  Aa ^ )  and  Cp(=C!  +  AC|),  is  intimately  related  to 
the  success  or  failure  of  the  image  restoration  process.  The  remainder 
of  this  section  is  devoted  to  demonstrating  the  inportance  of  these 
pole  loci  and  their  effects  on  the  restored  image  quality.  In  this 
chapter,  a  series  of  computer  processed  restored  images  is  generated  in 
order  to  determine  whether  the  presence  of  poles  does  indeed  influence 
the  allowable  ranges  for  the  parameters  Aaj  and  ac,.  Each  restored 
image  is  obtained  by  calculating  the  discrete  RIS  (Eq.  (59))  and  then 
applying  the  inverse  FFT  algorithm  to  this  function. 

Because  only  those  cases  where  aj  =  Cj  and  Aaj  =  ac}  are  being 
considered,  all  the  results  derived  throughout  the  remainder  of  this 
chapter  for  a]  and  Aai  apply  identically  for  C!  and  ACj. 

The  series  of  computer  processed  images  was  divided  into  three 
unique  cases.  For  each  case,  the  value  of  the  parameter  a,,  which 
describes  the  gross  wavefront  distortion  at  the  edge  of  the  input 
aperture,  was  specified  so  as  to  introduce  the  loci  of  poles  at  a  known 
location  in  the  RIS  (see  Fig.  10).  For  each  selected  value  of  a,,  a 
restored  image  was  calculated  for  perturbations  Aa]  =  n (0.006a)  (where 
n  +  10, 1,2, •••41).  The  special  situation  in  which  Aaj  =  0.0  corres¬ 
ponds  to  that  of  the  ideally  restored  image.  It  occurs  only  when  the 
degraded  image  and  the  PSF  have  experienced  exactly  the  same  turbulence 
effects  characterized  by  the  particular  choice  of  a!.  The  argument  Z, 
which  is  taken  from  Eq.  (58),  is  graphed  for  the  three  cases  in  which 
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31  °:5*  or]  °'75*  (equivalently,  in  terms  of  the  polynomial 

coefficients  introduced  in  Eq.  (8b)  these  three  cases  are  given  by 


2tt 


Cy/X 


m 


y(0.25x) , 


2tt 


Cy^X 


‘m 


y(0.50x)  or 


2tt  .  . 

F  7  ^(0  •  75  a  ) 
Son 


u^dein^FinaSn  °l  1n  Fi9*J2*,  graph>  which  is  similar  to  those 
used  in  Fig.  10,  presents  the  locations  of  the  pole  loci  in  the  RIS 

s/nd  °  nial?nDT?6  fis  re p resent  the  1,nitial  locations 

Thiel  *  S=  °-0)1of.RIS  P°le  loci  for  Cases  II  and  III,  respectively. 
Those  initial  pole  locations  are  defined  by  the  intersections  of  the 
Cot  argument  with  the  lines  1  =  mr. 


that  extend??™  ‘  ^  d’aracteristi':  <"*r  the  processing  range 
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(ai  -  0 . 246 )  <  'ap <  (a:  +  0.246) 


Fig.  12.  Location  of  the  points  at  which  Cot  argument 


"°  l°'l  !?d  f"‘°  the 
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SU.J;  11  iZ'lL  ‘actoMs  "graphed  as  a 

when  da,  =  o7o'(aD  .  f )  Sr  a°te  that  on  each 

are  the  expected  results  for  thP  *nd  1s  a  mi'ni'mum.  These 

all  graphs  independent  offthe  Choice  "for  a  1raa9e  a"d  °ccur  on 

ttBK&rSft^*  ss&ife 

:^sr»T;;irH  £  «'£.kv;s- ;,s.-:r 

r  and  £ms  are  monotonical ly  approaS^^"131"1"9  ?uality  Parameters 
point  Aa,  =  0.11a  Thic  5  9  ^eir  reactive  bounds  at  the 

that,  in  the  remaining  portion  of  the'Veainn1^  reinr(:ol[,ces  the  assumption 
rested  image  cannot  be 
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-0.246a  <  Aa,  <  0.11a 


where  -0.246a  is  the  lower  processing  limit  for  aa,  (see  Eg.  (70)). 

(a,  ^O^th'w^a^s-o'li/Jh^V4?  Ta,S  that  for  Ca^  1 1 
fms  vary  erratically.  Hidden  within  thesl°fl,!rtn  f?rameters.Sr,  r,  and 
in  which  the  restoration  navamofo  .  . fluctuations  we  find  regions 

define  a  successfully  restored  imaae  °  within  the  bounds  thal 

regions  in  which  the  oualiwJTnS'  Ho"?ver’  ‘here  are  adjacent 

bounds,  indicating  that  image  restoration  //ufor  Se'se^alue^of  Aa,. 

erratic9behaviorlorresponds  exactly^!  th^h^  *4c  ,indicates  that  this 
of  the  first  pair  of  pole  loci  ^  ^  introductl’on  into  the  RIS 

the  presence  of  polesPin  the  RIS  is  dirOn^^?0!^  the  hyP°thesis  that 
quality.)  Furthermore  when  Aa  <  n  l™  }i  ?*]***  to  restored  image 
exist  in  the  RIS.  This  explains  the  Jmnth?  indlcates  that  Poles  do  not 
parameters  in  Figs.  14a.  iSb  and  Hc.TJe’L^Mrl^0^  is 
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Fig.  13c.  ems  versus  Aaj  for  Case  I. 


versus  Aa5  for  Case  II. 


rersus  /.ai  for  Case  II 
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Fig.  14c.  ems  versus  Aaj  for  Case  II. 
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the  parameter  Sr  (see  Fig.  14a)  that  establishes  the  restoration  limit 
which  occurs  at  Aa1  =  -0.036 A. 

Thus  for  Case  II,  the  limits  on  Aax  are  a  series  of  isolated  small 
scale  regions  rather  than  large  scale  bounds  present  in  Case  I. 

Now  consider  Figs.  15a,  15b  and  15c.  These  graphs  indicate  that 
over  the  entire  range  of  Aaj  in  Case  III  (a.  =  0.75x),  there  are  only 
small  scale  regions  in  which  the  restoration  parameters  remain  within 
the  bounds  defining  a  successfully  restored  image.  Again,  these 
regions  are  defined  by  the  points  where  the  parameter  Sr  exceeds  its 
lower  and  upper  bounds.  However,  the  primary  difference  between  the 
small  scale  regions  in  Cases  II  and  III  is  their  width,  which  is  notice¬ 
ably  narrower  in  Case  II  (when  Aa^  >  0).  Another  basic  difference 
between  the  two  cases  is  the  initial  location  of  the  poles  within  the 
RIS  as  seen  in  Fig.  12.  In  Case  III,  the  poles  are  initially  (i.e., 
when  Aax  =  0.0)  further  removed  from  the  point  |KX|  =  Km. 

The  important  consequence  of  the  Case  II  and  Case  III  studies  is 
that  the  presence  of  RIS  pole  loci  imposes  limits  on  Aa,,  which  in  turn 
define  a  region  of  successful  image  restoration.  These  limits  are  a 
series  of  isolated  small  scale  regions  along  the  Aai  axis  in  Figs.  13a 
and  14a. 


Next,  an  expression  is  derived  that  relates  the  presence  of  RIS 
pole  loci  to  the  magnitude  of  Sr.  The  restored  Strehl  ratio  (Sr)  is 
easily  calculated  if  we  remember  that  Sr  is  the  spatial  domain  analog  of 
Wie  D.C.  value  in  the  frequency  domain  (see  Chapter  III,  section  E).  That 
is,  Sr  is  the  average  of  the  frequency  components  as  defined  in  Eq.(72a). 

63  63 

(72a)  Sr  =  £  R(p  AKx,q  AKy,  als  0,  Cj,  Aa^  0,  acj)  . 


Substituting  the  discrete  RIS  (obtained  by  substituting  Eqs .  (34a)  and 
(34d)  into  Eq.  (58)  with  a:  =  cx  and  Aa ^  =  acx)  into  Eq.  (72a)  we  obtain 


(72b)  Sr  =  4 


31 

I 

p=0 


Aa. 


1  + 


|(27r/Kln2)(2Knr|p  AKX |)  Aajp  aK> 


Cot  (2rr/Km2)  (2Km-|p  AKx  |  )ap  p  AKx 


The  terms  under  the  summation  sign  in  Eq.  (72b)  are  graphed  for  Cases 
I,  II  and  III  in  Figs.  16a,  16b  and  16c,  respectively.  These  graphs, 
which  are  the  envelopes  obtained  by  connecting  the  discrete  frequency 
components,  clearly  indicate  the  presence  of  poles  suggested  in  Fig.  12. 
Furthermore,  the  graphs  suggest  an  important  question.  How  is  the  value 
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rsus 


Fig.  15c. 


versus  Aaj  for  Case  III. 


of  Sr,  which  is  calculated  using  Eq.  (72b),  affected  by  whether  or  not 
a  pole  exists  exactly  at  a  discrete  frequency  component?  In  order  to 
answer  this  question,  a  relationship  is  derived  to  relate  the  magnitude 
o  Sr  to  the  possibility  that  a  discrete  RIS  frequency  component  was 
evaluated  at  a  pole. 


The  Cot  argument  Z(Kx,Ky)  taken  from  Eq.  (72b)  is  graphed  in 
Fig.  17a  for  Case  II  (aj  =  0.50x).  A  more  detailed  graph  of  this 
argument  in  the  neighborhood  of  the  point  Kx  =  Km  is  presented  in 
Fig.  17b.  (Fig.  17b  is  exaggerated  for  clarity.) 

Examining  Fig.  17b,  we  see  that  when  Aaj  =  0.0,  the  Cot  argument 
L  equals  the  value  ir  at  a  point  half  way  between  the  discrete  frequency 
components  15  aKx  and  16  aKx.  Thus,  since  the  ordinate  of  the  discrete 
argument  never  actually  equals  the  value  ir,  the  discrete  RIS  will  not 
contain  any  true  poles.  Hence  the  parameter  Sr,  which  is  calculated 
from  Eq.  (72b),  will  be  finite.  However,  referring  to  Fiq.  17b,  for 
the  curve  defined  by  Aaj  =  Aai',  the  Cot  argument  Z  does  equal  the 
value  y  at  15  aKx  and  16  aKx  (Km  ±  %aKx)  and  hence  RIS  poles  do  exist, 
or  this  situation,  the  summation  over  all  frequency  components  will 
obviously  exceed  the  limits  imposed  on  Sr.  When  the  value  of  an  aqain 
increases  to  the  value  aj  +  Aa^'  and  then  again  to  a,  +  Aa,"1,  the 
Cot  argument  behavior  is  analogous  to  that  of  Aa,  =  0  and  Aa,  =  Aa  1 
respectively.  That  is,  poles  do  not  exist  in  the  RIS  for  Aa!  =  Aa,1'' 

^hey  ?o  exist  Tor  Aaj  =  Aaj"1.  Thus,  based  on  this  analysis  and 
the  behavior  of  Sr  exhibited  in  Figs.  13a,  14a  and  15a,  it  seems 
natural  to  assume  that  small  scale  regions,  in  which  scccessful  image 
restoration  occurs,  will  exist  for  values  of  Aa!  in  ranges  similar  to 

fl.  -  ^r3)  Aai  .  (see  Fig.  17b).  That  is,  successful  restoration 
will  exist  in  regions  which  are  bounded  at  each  end  by  Sr  greatly  ex- 
ceeding  its  established  limits.  This  assumption  is  verified  in  Fiq  18 
which  is  a  more  detailed  graph  of  Sr  for  Case  II  (a.  =  0.50x)  where’ 

Aai  *  u  Irl  YVS  graph  many  sma11  sca1e  regions  are  apparent  and  each 
one  is  bounded  by  Sr  greatly  exceeding  the  limits  Sr  =  2.0  or  Sr  =  0.50. 

The  data  for  Fig.  18  was  obtained  by  using  the  same  computer  pro¬ 
cessing  techniques  employed  in  Cases  I,  II  and  III  but  with  finer 
resolution  in  the  parameter  Aaj  (i.e.,  Aaj  =  nO.OOl  ,  n  =  0,1,---  ). 

Next,  we  show  that  indeed  the  small  scale  regions  are  bounded  by 
those  values  of  Aa!  which  result  in  the  Cot  argument  Z  equalling  the 
value  tt  at  one  or  more  discrete  frequencies  causing  Sr  to  exceed  the 
limits  defining  a  successfully  restored  image.  In  order  to  prove  this 
fact,  an  equation  is  derived  that  specifies  the  magnitude  of  Aa, 
required  to  shift  a  qiven  ordinate  value  U  in  this  case)  of  the  Cot 
"n"mn  "  m/2  m  integer)  units  along  the  Kx  axis  (see 


argument  p  aKv  (p 
Fig.  17b). 


For  the  first  step  of  this  derivation  we  must  recall  that  the 
continuous  Cot  argument  (from  Eq.  (58))  is: 
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(73)  Z  =  (2,/y)  dp(  V  +  2 Kg.Kj,) 


Rewriting  Eq.  (73)  in  terms  of  the  discrete  frequency  components 
Ky  =  n  aKx  and  K,^  =  Ibh  aKx  (where  n  is  an  integer)  we  obtain: 

(74)  Z  =  (2tt/240 . 25 )  ap(-n2  +  31n) 

Eq-  (74)  is  next  expanded  about  the  point  Kx  =  Km  by  letting  [ pAKx 
(p  -  m/2,  m  -  integer)  measure  distances  along  the  abscissa  with  the 
point  Kx  -  Km  considered  as  the  origin.  Thus  the  Cot  argument  re¬ 
written  as  a  function  of  the  variable  p  is: 


(75)  Z  -  (2tt/240 .25 )  ap(240.25  -  p2). 


Finally,  substituting  aD  =  a,  +  Aa.  and  solving  for  Aa, ,  the  desired 
equation  is  obtained:  1 


(76)  Aa,  =  (Z/2T)(.?5° )  .  ,  . 

'240 .25  -  p2/ 

Eq.  (76)  specifies  the  magnitude  of  Aa,  required  to  shift  the  Cot 
argument  Z,  p  =  m/2  (m  integer)  units  along  the  Kx  axis. 

Applying  Eq.  (76)  to  Case  II,  in  which  a,  =  0.50A  and  Z  =  T  when 
A3j  =  0.0  (see  Fig.  12),  we  have:  1 

(77)  Aa,  =  0.5o( _ si _ )  . 

\240.25  -  p2/ 

Table  2  lists  the  values  of  Aa!,  calculated  from  Eq.  (77),  for  p  =  m/2 

■  =u-  u  ”’17  "  T^e  values  P  =  ni/2  (m  odd)  correspond  to  those  cases 
in  which  the  Cot  argument  exactly  equals  the  value  *  at  one  or  more 
discrete  frequency  (see,  for  example,  the  Aaj  =  Aai*  curve  in  Fig.  17b) 
and  poor  image  restoration  is  obtained.  The  values  p  =  m/2  (m  even) 
correspond  to  those  cases  in  which  the  value  tt  occurs  exactly  half  way 
between  two  adjacent  discrete  frequencies  (see,  for  example,  the 
Aai  =  Aai  curve  in  Fig.  17b)  and  good  image  restoration  is  obtained. 

Next,  a  series  of  computer  processed  images  was  generated  usinq 
the  same  processing  technique  described  earlier  for  Cases  I,  II  and  III- 
however,  the  values  for  Aax  were  taken  from  Table  2.  The  values  of  Sr 
calculated  for  each  restored  image  are  graphed  in  Fig.  19a  and  19b. 

Fig.  19a  is  a  superposition  of  the  data  obtained  from  Table  2  over  the 
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Fig.  19a.  Sr  evaluated  at  points  listed  in  Table  2 
and  superimposed  on  Case  II  Sr  versus 
A* i  graph. 
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Sr  evaluated  at  points  listed  in  Table  2 
and  superimposed  on  Fig.  18. 
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TABLE  2 


a  j  =  0.50A 


-  1 

p 

Aa  -  0.0005a 

=  1 

0.0021A 

~  1*2 

0.0047a 

=  2 

0.0084A 

-  2-k 

0.0133a 

=  3 

0.0194a 

—  3*2 

0.026 8 A 

=  4 

0  .0356A 

=  4*2 

0.0460A 

=  5 

0.0580a 

=  5k 

0 .0720 \ 

-  6 

0 .0881A 

=  6*2 

0.1067 A 

=  7 

0 . 1281 \ 

=  7*2 

0 . 1528A 

=  8 

0 . 1815  A 

-  8*2 

0.2150a 

Case  II  results  while  Fig.  19b  is  a  superposition  of  this  data  over  the 

D  -  JhSe  ’  d3ta  (i1"  =  n°-°01  data>'  Note  that  “hen 

L\„  „  (m  ?dd>  the  occurrence  of  the  Cot  „  situation  does  indeed  cause 

of  Sr  L  'ts  upper  bound.  When  p  =  m/2  (ra  even),  the  values 

thoS:  anr  flnite,  in  fact,  in  all  but  one  case,  shown  in  Fig.  19a, 

fully  res  tore  Vimagef1"  accePtab1e  bouads  for  Sr  defining  a  success- 

Finally,  we  derive  an  equation  that  expresses  how  rapidly  the  poles 
(i.e.,  points  corresponding  to  p  =  m/2  where  m  is  odd  in  the  above 
analysis)  migrate  through  the  RIS  as  a  function  of  Aa,.  This  point  is 
expressed  by  taking  the  derivative  of  the  Cot  argument  Z  with  respect 
to  Aaj .  The  result  of  this  calculation  is: 


(78) 


dl 


3  Aa  i 


=  (2r/V)(2Km  - 


m 


Kx!)K> 


wMch  is  graphed  in  Fig.  20  as  a  function  of  Kx.  Fig.  20  demonstrates 
that  for  a  given  perturbation  Aa1}  the  poles  mf grate  from  one  discrete 

c?!ipo,?enJ  t0  t!?e  next  adjacent  component  most  rapidly  (i.e., 
o  the  smallest  change  in  Aai)when  they  occur  near  the  point  K  =  K 

/,1S  fnCcnC^rTSP°nds  exactly  with  the  behavior  seen  in  Case  II  X 
ia!  T.0:50^  wh^  rapid  fluctuations  occur  in  Sr  as  Aa,  increase  in 

oblerve^in'case  in  ^  “c  ^  Fig-  20  a1s°  expla1'n^  the  behavior 
°a  l065  not  Ya^  as  rapidly  as  a  function 


nr 


cn„  r„_  TTT  /r-.:.:  i  ap i iy  ds  a  Tuncnon 

Case  III  (Fig.  14a)  when  Aa2  is  positive,  the  poles  migrate 
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Fig.  20.  Graph  of  versus  KY. 

x 

farther  away  from  the  point  Kx  =  1^  and  the  small  scale  regions  have 
larger  separations;  however,  when  Aa,  is  negative,  the  opposite  effect 
is  observed. 

We  can  now  generalize  from  these  results.  At  best  we  can  expect 
only  small  scale  regions  in  which  the  imaging  processing  technique  will 
produce  a  successfully  restored  image.  In  fact,  from  the  data  presented 
in  Figs.  13a,  14a  and  15a  we  can  calculate  the  percentage  of  the  Aaj 
range  over  which  successful  restoration  results  occur.  For  the  three 
cases,  these  percentages  are: 

CASE  I  85.5% 

CASE  II  14.0% 

CASE  III  71.0%  . 

The  presence  of  these  small  scale  regions  will  obviously  create  great 
difficulties  in  the  analysis  of  the  restored  images. 

The  fact  that  a  particular  image  was  not  successfully  restored 
could  be  the  result  of  one  of  two  factors.  First,  the  object  and  point 
source  could  actually  have  been  located  in  different  isoplanatic  patches. 
By  definition  of  an  isoplanatic  patch  (see  Chapter  I )  we  would  expect 
poor  quality  restoration  in  this  case.  Second,  the  atmospheric  effects 
which  degraded  the  image  might  have  created  a  situation  in  which  the 
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Cot  argument  in  the  RIS  assumed  the  value  v  at  one  or  more  discrete 
frequency  components.  As  we  have  seen,  this  event  will  grossly  perturb 
the  RIS  so  that  when  it  is  inverse  Fourier  transformed,  the  restored 
image  does  not  meet  the  restoration  criterion  established  in  Chapter 
III,  section  F. 

C.  Summary 

This  chapter  has  been  devoted  to  the  analysis  of  three  cases  of 
restored  images.  The  three  cases  were  defined  by  the  values  assumed 
for  the  wavefront  parameter  aWq  =  a^.  Each  case  was  analyzed  for 
a  series  of  perturbations  Aa1(Ac1  =Aa1)  in  the  initial  value  for  a:. 
The  restored  image  quality  parameters,  restored  Strehl  ratio  (Sr), 
integral  scale  (r),  and  mean  squared  error  (e^)  were  calculated  for 
each  restored  image,  and  then  graphed  as  a  function  of  Aax.  From  these 
graphs,  the  regions  of  successful  image  restoration  are  defined  between 
the  points  at  which  Sr  exceeds  its  upper  and  lower  bounds.  The  exceed¬ 
ingly  large  values  for  Sr  were  then  related  to  whether  or  not  a  pole 
was  evaluated  exactly  at  a  discrete  RIS  frequency  component. 

Finally,  it  was  demonstrated  that  the  number  of  small  scale  regions 
for  a  given  restoration  range  (i.e.,  range  of  ap)  is  dependent  on  the 
initial  location  of  the  poles.  For  poles  initially  located  near  the 
point  Kx  =  Km,  the  poles  migrate  rapidly  from  one  discrete  RIS  frequency 
component  to  the  next  adjacent  frequency  component.  Thus,  the  regions 
of  successful  image  restoration  as  a  function  of  Aa:  are  very  narrow 
(as  small  as  0.0005A).  For  poles  initially  located  at  a  distance  from 
the  point  Kx  =  Km,  larger  variations  in  Aaj  are  required  to  shift  a 
pole  from  one  discrete  frequency  to  the  next  adjacent  component.  Thus, 
the  regions  of  successful  image  restoration  encompass  a  larger  range 
of  Aax. 

The  important  point  established  in  this  chapter  is  that  image 
restoration  is  a  much  more  complex  process  than  dividing  a  degraded 
image  spectrum  by  a  MTF  and  then  inverse  Fourier  transforming  to  obtain 
the  restored  image.  It  has  been  shown  that  the  poles  introduced  by 
this  division  operation  limit  image  restoration  to  a  series  of  small 
regions  defining  perturbations  in  the  reference  point  wavefront  that 
are  not  present  in  the  point  image  wavefront. 
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CHAPTER  V 

REVIEW  AND  CONCLUSIONS 


A.  Review  and  Conclusions 

This  study  has  been  devoted  to  the  task  of  defining  the  limits  for 
which  the  Fourier  transform  method  of  image  restoration  may  be  used 
to  successfully  restore  turbulence  degraded  images  obtained  from  a 
ground-based  telescope. 

In  order  to  define  the  limits  over  which  successful  image  restora¬ 
tion  could  occur,  a  turbulence  degraded  imaging  system  was  modeled,  with 
the  assumption  that  the  degrading  effects  could  be  expressed  by  expand¬ 
ing  the  E-field  phase  distribution  across  the  telescope  input  aperture, 
in  terms  of  a  general  second  order  quadratic  equation.  This  expansion’ 
was  used  to  define  the  degraded  E  field  for  a  point  object  and  the 
degraded  E  field  from  a  reference  point,  which  wa^  assumed  to  be 
located  in  the  transmitter  plane  with  the  point  object.  In  order  to 
model  various  reference  point  -  object  point  separation,  the  reference 
point  E-field  phase  distribution  included  a  perturbation  term  in  the 
quadratic  phase  coefficients. 


The  actual  restored  image  analysis  discussed  in  this  study  was 
based  on  a  simplification  of  the  general  quadratic  phase  distribution 
discussed  above.  The  specific  case  analyzed  included  only  the  spherical 
wavefront  distortion  in  the  E-field  phase  distributions.  With  the  input 
aperture  considered  as  a  spatial  frequency  plane  (see  Chapter  n-B)  the 
object  point  and  reference  point  E- field  phase  distributions  for  the 
spherical  wavefront  distortion  case  are,  respectively: 


(79a) 

(79b) 


,  .  .  J^A)2  Mx2  +  c l I'y  ) 

Ea(Kx,y  «  e  1  y 

(  j(di/k)2  {(ai  +  Aal)K2  +  (Cl  +  ACl)Ky2  } 

a  x’  y  e 


The  parameters  a1  and  cl  represent  gross  wavefront  distortions  while 
the  parameters  Aa}  and  ac:  define  the  phase  perturbations  resulting 
from  the  object  point  and  reference  point  waves  propagating  through 
different  turbulence  conditions.  (See  Chapter  II,  section  B  and  section 
C  for  details.) 

Eqs.  (79a)  and  (79b)  were  used  as  inputs  for  calculating  the 
degraded  image  spectrum,  and  the  modulation  transfer  function  (MTF) , 
respectively.  Then,  using  linear  system  theory,  the  restored  image  was 
obtained  by  calculating  the  inverse  Fourier  transform  (The  Fast  Fourier 
iransform  (FFT)  was  used  for  computer  processing  (see  Appendix  C)  of  the 
degraded  image  spectrum  divided  by  the  MTF  or  system  transfer  function  ) 
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The  restored  image  spectrum  (RIS)  obtained  by  dividing  the  de¬ 
graded  image  spectrum  by  the  MTF  was  subjected  to  a  complete  analysis. 

It  was  shown  that  the  presence  of  poles  (points  at  which  the  mathe¬ 
matical  function  takes  on  an  infinitely  large  value)  in  the  RIS 
definitely  affected  the  restored  image  quality,  fquations  were 
derived  that  related  the  number  and  position,  within  the  RIS,  of  these 
poles  to  the  parameters  a]  +  /\al  and  Cj  +  ACj . 

The  restored  image  quality  was  defined  in  terms  of  three  resto¬ 
ration  parameters:  the  restored  Strehl  ratio  (Sr)  which  measures  the 
peak  image  intensity,  the  integral  scale  (r)  which  measures  image 
width,  and  the  mean  squared  error  (rms)  which  are  discussed  in  detail 
in  Chapter  III,  section  E.  Limits  were  placed  on  these  parameters 
which  defined  a  successfully  restored  image. 

A  series  of  degraded  images,  generated  using  the  model  discussed 
above,  was  restored  using  a  digital  computer  to  implement  the  inverse 
Fourier  transform  restoration  method.  The  degraded  images  were  divided 
into  three  classes  defined  by  the  value  assumed  for  the  parameter 
ai(=  ^i)-  For  each  such  class,  a  series  of  restorations  were  performed 
for  perturbations  Aa1  =  n(0.006A)  (n  =  ±0,1, 2, •••41). 

The  results  of  the  analysis  of  these  restored  images,  in  terms 
of  the  restoration  parameters,  indicated  that  the  presence  of  poles  in 
the  RIS  did  indeed  cause  the  restored  image  to  fail  to  meet  tne  re¬ 
stored  image  quality  criterion.  It  was  determined  that  the  quality 
parameter  Sr  was  the  most  sensitive  indicator  of  restored  image  quality, 
and  hence  this  parameter  was  used  to  define  the  ranges  of  successful 
image  restoration. 

The  conclusion  drawn  from  this  study  was  that  a  discrete  RIS 
frequency  component  evaluated  precisely  at  a  pole  was  the  major  factor 
that  caused  Sr  to  greatly  exceed  its  bounds  indicating  that  the  restored 
image  failed  to  meet  the  quality  criterion.  As  these  poles  migrated 
from  one  discrete  spatial  frequency  component  to  the  next  adjacent 
discrete  component,  a  region  of  successful  image  restoration  was 
bounded.  These  small  scale  regions  were  found  to  be  narrower  in  width 
(the  width  expressed  in  terms  of  Aa  )  when  the  poles  originated  near 
the  spatial  frequency  spectrum  midpoint  Kx  =  K  Poles  initially 
located  some  distance  from  this  point  Kx  =  Km  were  found  to  generate 
much  wider  small  scale  regions  of  successful  image  restoration. 

For  the  three  cases  studied,  it  was  shown  that  the  sum  of  these 
small  scale  regions  ranged  from  14%  to  85.5%  of  the  processing  range 
that  was  considered.  Obviously,  the  14.0%  case  indicates  a  very  severe 
limit  imposed  on  the  Fourier  transform  image  restoration  method. 
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The  results  obtained  from  this  study  support  the  following  six 
conclusions: 

1.  The  major  problem  in  the  restoration  of  images  is  the 
presence  of  poles  in  the  restored  image  spectrum  (RIS). 

2.  The  restored  Strehl  ratio  (Sr)  is  the  most  sensitive 

of  the  three  restoration  parameters  (integral  scale,  r, 
and  mean  squared  error,  ,  are  the  other  two  restora¬ 
tion  parameters)  used  to  measure  restored  image  quality. 

3.  The  ability  to  successfully  restore  a  degraded  image 
does  not  monotonically  decrease  with  increasing  wave- 
front  distortion  but,  rather,  is  oscillatory  in  nature. 

4.  Given  a  gross  wavefront  distortion,  there  may  be  many 
small  ranges  of  wavefront  perturbations  in  which 
successful  image  restoration  occurs. 

5.  To  assure  restoration  in  the  al  =  0.50x  worst  case 

considered  here,  Aa.  must  be  less  than  or  equal  to 
0.0005X.  1 

6.  A  orobabilistic  measure  of  successful  image  restora¬ 
tion  as  a  function  of  gross  wavefront  distortion 
appears  to  be  a  useful  discriptor. 

B.  Suggestions  for  Further  Study 

In  this  section,  we  propose  suggestions  for  the  possible  extension 
of  the  work  described  in  this  study. 

The  next  step  beyond  the  work  covered  here  might  be  that  of  em¬ 
ploying  a  smoothing  function  such  as  a  two-dimensional  Gaussian  or 
triangular  function  to  smooth  out  the  effects  of  poles  present  in  the 
restored  image  spectrum  (RIS)  and  hence  improve  the  restored  image 
quality.  Further,  this  work  suggests  the  need  for  investigating  the 
validity  of  the  discrete  representation  for  the  quotient  of  two 
functions  as  in  calculating  the  RIS.  The  basic  question  to  be 
answered  is:  do  poles  in  the  RIS  result  in  effects  such  as  a 
spatial  domain  analog  to  aliasing? 

The  other  important  extension  would  be  that  of  including  the  cross 
term  parameters  bj  and  Abj  in  the  model  for  the  degraded  E-field  wave- 
front.  We  would  then  be  interested  in  the  location  of  the  new  RIS  pole 
loci  and  the  effects  of  these  poles  on  the  size  of  the  regions  in  which 
successful  image  restoration  is  obtained. 
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APPENDIX  A 

THE  APERTURE  PLANE  CORRELATION  FUNCTION 


field’aperture'distHbutio^fEq^tM^J^itMtself0^0^^0^^6^  I’t 

plane  correlation  function  C(KX,K  )  the  dperture 

a,  sar=  s  S'fs.i'js.arr^B 

IJI.JIJJ.  ,»,.l  Fourier  Mo,,.™  S',,*,™ 


(20) 


I(x2,y2)-  (1/2tt)2  ffu(Ky  _  k 


VX  >Ky 


W»(-k 


X-'ky) 


Ea(K> 


ky  >  K. 


U  dkY  dk, 


U9)  W(Kx,Ky)  =  ) 

,0 


K 


"Km  ^  Kx<- 

*m*  K;<  K 


m 

m 


otherwise 


we  have  llll'nln%\  a"d  "lultiplyin9  E<>-  <20>  by  (2’>2, 


where  C(Kx,Ky)  is  the  aperture  plane  correlation  function. 


We  now  are  interested  in  evaluating  Eq.  (21)  for  the  casP  whprp 
the  input  field  distribution  is  given  by  Eq.  (14)  as: 


(14) 


Ea(Kx»Ky) 


the 


Substituting  Eq.  (14)  into  Eq.  (21),  we 
actual  integral  that  we  must  evaluate. 


obtain  Eq.  (23)  which  is 
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(23) 


j(di/k)2(aKx2  +  bKxKv  +  cK v2) 

C(Kx,Ky)  =  e  x  y  y 

ff  J  ( d-s/k ) 2 

JJ  W(Kx+kx,Ky+ky)  W*(kx,ky)  e 

-  oo 

[(2aKx  +  bKy)kx  +  (2cKy  +  bKx)ky] 

^  J  uKy 


Since  the  aperture  function  restricts  the  input  field  at  ±Km, 
we  can  determine  the  limits  of  integration  by  considering  the  correla¬ 
tion  of  the  two  aperture  functions.  Fig.  Al,  below,  shows  the  two 
aperture  functions  involved  in  the  Kx  correlation. 


Fig.  Al. 


The  complex  conjugate  notation  {  *  )  may  be  dropped  from  W  (kx,k  )  since 
the  aperture  function  is  a  real  function.  Using  Fig.  Al  to  define  the 
limits  of  integration,  we  find  that  when  Kx  >  o)  the  limits  on  kx  are 

(Ala)  -Km  <  kx  <  -Kx  +  Km 

and  when  Kx<0,  the  limits  jn  it  are 

(Alb)  -Kx-Km<  kx<  Km 
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Since  the  input  aperture  is 
correlation.  That  is,  the  ky 


square,  the  same  results  apply  for  the 
integration  limits  are  defined 


(A2a)  for  Ky  >  0 

(A2b)  for  Ky  ^  0 


-Km  <  ky  <  -Ky  +  and 

-Ky-Km<  Ky  <  Km 


i-ntegrn1^Eq-E|ch3^"egbrat^r]stth,»„aS/h^PJr0!IUCt  0f  a  k*  and  2  ky 

by  Eqs.  (Al)  and  (A2).  9Note  that  the  an1Vfded  lnt0  tW0  parts  defl'ned 

since  its  effect  s  now  reflected  S  the  fu,Jctl0n  Can  be  dr<W 
we  have:  rejected  in  the  limits  of  integration.  Thus 


(A3) 


*  y 


—  ■ »» 


C(Kx,Ky)  =  e 

f-Kx+Km  K 

j(dj/k)  2(2aKx+bKy  )kx  51 

6  dkx  +  I  e 

"Kx_Km 


m 


j(di/k)2(2aKx+bKy)kx 

dk. 


"Ky+Km  K 

i  J  ej(di/k)2(2cKytbKx)kyd|(  +  ? 

y  dkJ 


For  the  first  term,  when  Kx  >  0,  we  have: 


(A4) 


~/X+^m 

f  eJ(d1./k)2(2aKx+bKy)kx 


-K 


dk„  = 


j(dj/k)2(2aK  +bK  )  Kx+Km 
e  *  y 


"  j(dj/k)2(2aKx+bKy) 

j(d1  /k )  2  ( 2a  Kx+bKy )  (-Kx+Km)  j(d1-/k)2(2aKx+bKv)K. 


-K 


m 


y  /INm 


j(dj/k)2(2aKx+bK  ) 
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j  (df/k  )2!j(2aKx+t. Ky  )KX 


d,-  /k  )2  (2  a  Kx+b  Ky )  2  j 


j  (di/k)2(2aKx+bK Jd^yy 


j(  d,.  /k ) 2  (2a  Kx+bKy )  ( K,,,-  %KX ) 


Eq-  (A4)  Can  be  ' 


in  terms  of  the  Sin 


,+K, 


m 


(A6) 


j(d1-/k)2(2aKx+bKy)k) 


dkx  = 


^di/k)^(2aKx+bKy)Kx  Sin(di/k)2(2aKx+bKy)(Km-J5Kx; 


%(di/k)2(2aKx+bKy) 


For  the  second  Kx  term  in  Eq.  (A3)  we  have  for  Kv  <  0* 


*Sn 

(A7,  J 

-Kx-Km 


j(2aKx+bKy)k. 


e  *„  -  - 


j(2aKx+bKy)kx 


j(2aKx+bKy) 


■Kx-Km 


j(d1/k)2(2aKx+bKy)Km  j(d1-/k)2(2aKx+bKy)(Km+Kx) 


j(d1-/k)2(2aKx+bKy) 


-j(d17k)2J5(2aKx+bKy)Kx 


^(dj/k)2 (2aKx+bKy)2j 


_J!  ^  j(di/k)2(2aKx+bKy)(Km+}5Kx) 

-j(d1./k)2(2aKx+bKy)(Kftl+}5Kx) 
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Again,  this  last  equality  in  Eq. 
Sin  function  as: 


( A7)  can  be  written  in  terms  of  a 


( A8 ; 


-K  +K 

/  j(di/k)2(2aK  +bK  )k 

e  *  y  »  dkx  = 


-j (d-/k)2^(2aK  +bK  )K 
e  1  x  y  x 


[Sin(d17k)2(2aKx+bKy)(Km+isKx]r| 


L  Js(dj/k)2(2aKx+bKy) 


Thus  for  the  Kx  correlation  we  have  from  Eqs.  (A6)  and  (A8)  that: 

K„ 


-Kx+Km  . 


(A9) 


j(d1-/k)2(2aKx+bKy)kx 
e  dkx  + 


'm 


j(di/k)2(2aKx+bK  )kx 


dk. 


-Km 


"^x"Km 


j(di/k)2J5(2aKx+bKy)Kx|pin(d1-/k)2(2aKx+bKy)Js(2Km-Kx) 


J5(di/k)2(2aKx+bK  ) 


^  Sin(di/k)2(2aKx-t-bKy)4(2Km+Kx)> 
!5(dj/k)2(2aKx+bKy ) 


We  next  multiply  the  first  Sin  term  (defined  for  Ky 
^tn'Kx  and  second  Sin  term  (defined  for  Kx  < 
(2KmtKx)  whl‘ch  casts  these  terms  in  the  form  of  the 
or  Sine  function.  Thus  the  right-hand  side  of  Eq. 


>  0)  by  (2K  -Kx)/ 
0)  by  (2Km<)/X; 
familiar  Sin  x/x 
(A9)  becomes: 
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(A10) 


j  (di/k)2J5(2aKx+lKy)Kx 


Sin(d-i/k)2(2aKx+bKy)^(2Km-Kx) 

(di/k)2(2aKx+bKy)J5(2Km-Kx) 


+  (2Km+Kx) 


Sin(di/k)2(2aKx-(-bKy)is(2KmrKx ) 

(di/k)2{2aKx+bKy)Js(2Km+Kx) 


Eq.  (AlO)  may  be  writion  in  a  more  compact  form  by  employinq  |KVI  ard 
reme^HHg  that  the  first  tens  in  Eq.  (AlO)  is  (tefin^d  for  K  >  0 
and  the  second  term  is  defined  for  Kx  <0.  x  “ 


(All) 


j(d-7k)2Js(2aKx+bKjKx 
e  {2Km-|Kx|) 


Sin  (di/k)2(2aKx+bKy)^(2Km-|Kx|) 
(di/k)2(2aKx+bKy)55(2Km-|Kx|) 


substituting-11  ^  S°1ution  to  the  oecond  bucketed  term  in  Eq.  (A3)  by 


Kx  for  Ku 
Ky  for  Kj 
c  for  a 

into  Eq.  (All)  to  obtain: 

(A12)  =  ej(di/k)^2cKy+bKx)Ky  (2V||<  [}  Sin(dj/l^(2cKy+bKx)^(2Kni- |Ky | ) 

y  (di /k )  2(2cKy+bKx)Js(2Km- 1  Ky  | ) 


Substituting  Eq.  (All)  and  (A12)  into  Eq.  (A3),  we  see  that  the 
multiplicative  phase  term  is  cancelled  leaving  us  with: 


(A13) 


C(Kx,Ky)  =  (2Km-|Kx|) 


Sin(di/k)2(2aKx^bKy)^(2Km-|Kx|) 
(di/k)2(2aKx+bKy)Js(2Km-  |KX|) 


(  |K  d  $in(d^/k)2(2cKy+bKx)is(2Km-|Ky|) 

y  (di/k)2(2cKy+bKx)!5(2Km-|Ky|) 
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Eq.  (A13)  which  is  defined  for: 


(AM) 

-2Km  i  Ky  <-  »■ 

is  the  desired  result  and  is  used  in  Eq.  (24)  of  Chapter 


APPENDIX  B 

CALCULATION  OF  THE  RESTORED  IMAGE  SPECTRUM 


Td/u  In,  aPPendix,  we  calculate  the  restored  image  spectrum  Iris) 

spectrum  C?KVe|<  {n,Chapter  "•  Eq*  (29^  the  ^graded  fZe^  ^ 

modulation  t?ansfe^Ccti^nP(MTF)CH[r1Kt)0n  fUnCtion)  and  the  system 

'  '  x*  y  ’ 


(29)  IR(Kx,Ky)  = 


C(Kx,Ky) 


H ( Kx , Ky ) 


|2t)  „and  3lvfr..:^  -ChaPier 


c  /nA\  -  .  .  emu 

Eqs  (24)  and  (28),  respectively.  The  fbtlowing  derivation  "of"the“RIS 

:iVrl«r,veTK.onlLlh.e,  2“  >y.k:|K  *1)"  bM) 

/I)  is  ( 


and  (28).  The  second  Sc  tern muTi ed by  (2kS  fflW:  (H 
directly  from  the  following  deri™^ l? XV Sibsii  tStiorlsf 


(Bla) 

and 

(Bib) 

and 

(Blc) 


Ky  for  Kx 


cx  for  a: 


ACj  for  Aax 


in  theJRIS:tUt’"9  theSe  S'"C  tennS  1nt°  Eq‘  (29)’  we  have  as  °"e  tern, 


(B2) 


Sin{(d1-/k)2(2aKx+bKy)^(2Kn|- 1  Kx  I )} 
(df/k)2(2aKx+bKJis(2K  -|k  |) 


Sin{(d^/k)2f2(a+Aa)Ky+(b+ab)Kv]?s(2Km-)Kv 
(d17k)2t2(a+Aa)Kx+(b+4b)Ky]!s(2K  -|K„|) 


(is" i  SXdsre*^rtl!?Siiit4(Vk)a(2l?/fh*t&Adkf,2!3s-  (15?>  a?„d,, 

as  a  ratio  of  two  Sin9  funalo^it^TeJ  b/l%ioTiZ?^X  ind’ 


■Xy. 
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A 


Sin  |(27r/Km2 )  (2a  iKx+bI K^,)%( 2 Km- 1  K„  | ) | 

Sin{(2»/Km2)[2(a1+4a1)Kx+(b1+4b1)Ky]%(2Km-|Kx|)|, 

We  next  add  and  subtract: 

(B4)  aY  =  (2„/Km  )(2Aa,Kx  +  Ab1Ky)'J(2Km-|Kx|) 

whi'ch ‘we  STSTL f  the  Si"  ,UnCti°n  fn  the  '’Umerator  °f  E<l-  (B3) 
(B5)  Y  =  (2,/Km2)(2a,Kx  +  b1Ky)fst2Km-|Kx|>  . 

Thus,  Eq.  (B3)  becomes: 


2Aa1Kx  +  AbjKy  I  Sin(Y+AY-AY) 
2aiKx  +  bjKy  .  Sin(Y+AY) 


Applying  the  trigonometry  identity. 

Sin  A-B  =  Sin  A  Cos  B  -  Cos  A  Sin  B  , 
to  Eq .  (B5),  we  have: 

2AaiKX  +  AbtKy  1  Si n ( Y+ aY )  Cos  AY  -  Cos(Y+aY)  Sin(AY) 
2aiKx  +  biKy  J  Sin ( Y+a Y ) 


which  simplifies  to: 


+  2Aai*S<  +  Ab^/ 

2aiKx  +  bjKy 


(B7) 


Cos  AY  -  Sin  aY  Cot(Y+AY) 


Substituting  Eqs.  (B4)  and  (B5)  into  (B7),  we  have: 


(B8) 


r  2Aa1Kx  +  Ab,K  I 

1  +  - — - H 

L  2alKx  +  blKy  J 

Ccs  |(2ir/Km2)(2iaIKx  -  ib1Ky);s(2Knl-|Kx|)l  - 
Sl'n|(2Ti/K|n2){2aa1Kx  -  Ab,Ky)>s(2Km-|Kx|)} 

-  Cot-ffa^/l^z) [2(a1+AaI)Kx  +  (b1+ab1)Kyj!s(2Kn)-|Kx|)} 


As  stated  earlier,  the  second  term  in  the  RIS  is  obtained  by  making  the 
substitutions  indicated  in  Eqs .  (Bla),  (Bib)  and  (Blc)  into  Eq.  (B8)  to 
yield:  M  v  ' 


(B9 ) 


2  +  2AciKy  +  Ab!Kx 

2clKy  +  biKx 


[^{(^/VX^CiKy  -  4blKx)y2Km-|Ky|)}  - 

Sin{(2„/y)(2ac1Ky  -  ib,Kx)%(2Km- |Ky | )| 

•  Cot  |(27r/Km2)[2(a1+AC1)Ky  -  (b1+Ab1)KxJJs(2Km-|Ky|)J. 

The  desired  RIS,  R(KX,K ,a, jbj ,c, ,Aaj  ,Abj ,ACj) ,  in  Eq.  (31)  of  Chapter  n 
is  the  product  of  Eqs.  {B8J  and  (B9).  H 

(BIO)  R(Kx.Ky,a1,b1,c1.Aa1.4b1.Ac1)  =  [l  ♦  2_Aa»Kx  j  abiKy 

L  Ba^+AbjKy. 

|cos{(2./Kn2)(2Aa1Kx  -  4b,Ky)!5(2Knl-|Kx|))-Sin/(2„/Kn12) 
•(ZAa1Kx-Ab,Ky)%(2Kn,-|Kx|)j.Cot/(2»/Km2)[2(a1+Aa1)Kx-(b]+Ab1)Ky] 


2ac.K  +  Ab,KY 
1  +  — LI - * 

2ClKy  +  biKx  _ 


Cos|(27r/Km2)(2Ac1Ky-Ab1Kx) 
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APPENDIX  C 

THE  FAST  FOURIER  TRANSFORM  (FFT) 


In  this  appendix  we  discuss  pertinent  aspects  of  the  Fast  Fourier 
Transform  (FFT)  algorithm  and  the  Discrete  Fourier  Transform  (DFT)  from 
which  it  is  derived.  Attention  is  also  given  to  both  frequency  and 
spatial  domain  discrete  variables  and  their  relationships  with  the 
spatial  domain  sampling  rate. 


It  is  a  well  known  fact  that  for  a  frequency  band-limited  system, 
the  time  waveform  is  exactly  described  by  a  Discrete  Fourier  Series 
(DFS),  when  the  time  function  is  sampled  at  a  rate  greater  than  or 
equal  to  the  highest  frequency  component  in  the  pass  band.  This 
relationship  applies  equally  for  a  n-dimensional  system.  The  Discrete 
Fourier  Transform  (DFT)  applied  to  the  DFS  is  a  bi-directional  mapping 
operation  with  mathematical  properties  analogous  to  tnose  of  the 
Fourier  integral.  The  Fast  Fourier  Transform  (FFT)  is  a  method  for 
efficiently  computing  the  DFT  of  a  n-dimensional  series  of  discrete 
data  samples.  For  the  purpose  of  this  study,  we  restrict  our  attention 
to  the  two-dimensional  FFT.  The  FFT  takes  advantage  of  the  fact  that 
the  coefficients  for  the  DFT  may  be  calculated  by  using  a  parallel 
iterative  technique  rather  than  the  coefficient-by-coefficient  direct 
approach.  For  a  two-dimensional  array  of  dimension  N  by  N,  where 
N  =  2m  (m  integer),  direct  calculation  of  the  DFT  coefficients  requires 
N4  arithmetic  operations;  however,  the  FFT,  using  the  'terative  tech¬ 
nique,  requires  only  4N2log,N  =  4N2m  arithmetic  operations  [181.  When 
N  =  64,  as  in  this  study,  the  FFT  requires  only  4m/N2  =  24/4096  or 
0.6%  of  the  possible  N4  arithmetic  operations,  and  hence  computer  time, 
required  by  the  direct  method.  Not  only  does  the  FFT  reduce  the  re¬ 
quired  computer  time  by  99.4%,  but  it  reduces  the  computer  round  off 
errors  by  the  same  factor  [*18] . 


The  defining  equation  for  the  DFT  is: 


(Cl) 


IJ 


.  Nx  Ny  f  __  ♦  liljjfidl, 


R=1  S=1 


where  IgS  is  an  Nx  by  Ny  array  of  discrete  data  samples  and  Ijj,  the 
DFT  of  Inc,  is  an  Nx  by  Ny  array  of  discrete  frequency  components.  The 
inverse  DFT  is  given  by: 


(C2) 


lRS 


_  1 


Nx 

v 

L 


Vy  .=1  J: 


i  2tt  ( 


(I-1)(R-1)  (J-D(S-l) 


Nx 


IJ 


Ny 


) 


The. validity  of  this  transform  pair  can  be  easily  verified  by  substi¬ 
tuting  Ijj  from  Eq.  (Cl)  into  Eq.  (C2)  and  using  the  orthogonal 
relationship  for  exponentials,  which  is: 
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(C3) 


Nx 

l 

1  =  1  J=1 


Ny  _  s2tt( 


(I-1)(R-1)  (J-1)(S-1) 


Nx 


Ny 


( 1-1 ) (M-l )  (J-1)(N-1) 

>  J2"(  Mv"  ' +  ~~ N7 — > 


Nx 


(  Nx  Ny 


1o 


for  R  =  M 
S  =  N 


otherwise 


For  convenience,  the  DFT  transform  pair  is  represented  by  the  shorthand 
notation: 


(C4)  IRS  Ijj 

We  are  assuming  in  this  study  that  the  frequency  spectrum  is  band- 
limited  at  Kx  =  ±2Km  and  Ky  =  ±2Km  (see  Eq.  (21));  thus  the  sample 
spacing  Ax  =  Ay  in  the  spatial  plane  must  be  given  by: 


(C5)  ax  =  Ay  <  %(2ir/2Kra)  =  (u/2Km) 


We  further  assume  that  the  width,  W,  of  spatial  function  equals  an 
integral  number  of  samples;  that  is: 


(C6 ) 


W 


Nxir 

NxAx  =  NyAy  -  - 


Nyir 


Substituting  Nx  =  W/ Ax  and  Ny  =  W/Ay  into  Eq.  (Cl)  and  (C2)  and  grouping 
terms  in  the  exponents,  we  have: 


(C7 ) 


IJ 


Nx  Ny 

I  J  1 RS  e 
R=1  S=1 


(I-D 


(J-l) 


-j2,r("ir<R-1)Ax  + 


and 

(C8) 


I 


RS  "  N  N 


(I-D  (J-l) 

i  N*  (R-»4X  *  -jj-  (s-i)iy) 

•j  w  L  l  ^IJ  e 
xy  1=1  J=1 
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From  Eq.  (C7),  we  see  that  is  periodic  with  period  l/&x  along  the 
Kx  axis  and  1/Ay  along  the  Ky  axis;  that  is: 

<C9>  >U  -  W«*.J±n/4y  ("  1nte9er)  ‘ 

Alsb,  we  see  that  the  sample  spacing  in  the  frequency  plane  is: 


(CIO)  AKX  =  AKy  =  1/W  . 

From  Eq.  (C8),  we  find  that  Ir$  is  periodic  with  period  W  along  the  x 
and  y  axes ;  that  is : 


<C11>  JRS-WstnH  ("integer)  • 

The  sample  spacing  has  already  been  specified  as  Ax  =  Ay.  Further, 
note  that  since  the  spatial  domain  period  W  =  NxAx  =  NyAy  (from  Eq. 
(C6))  and  the  frequency  plane  periods  1/Ax  =  Nx/W  and  1/Ay  =^Ny/W 
(from  Eq.  (C6)  are  both  directly  proportional  to  Nx  and  Ny,  I j j  and 
1^5  are  also  periodic  with  periods  Nx  and  Ny.  Thus,  each  array  of 
Nx  by  Ny  number  represents  one  cycle  of  the  periodic  frequency  and 
spatial  domain  functions. 

These  results  are  summarized  in  Fig.  Cl  for  the  one-dimensional 
case  (i.e.,  Ij  <h  Ir),  which  simplifies  the  graphic  presentation. 
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